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Nonlinear  optical  snsceplibililies  provide  a  convenient  means  of  relating  macroscopic  optical  measurements  to  micros^oric  nioocis  1  he 
susceptibilities  .ire  most  uselul  when  the  radiation  field  and  the  material  degrees  ol  Ireedom  are  we.ikiv  coupled.  In  the  opposite  w.ise.  ilic  I'.vn.imies 
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1.  Introduction 

The  transient  grating  (TG)  technique  is  one  of  the  most  powerful  means  for  the  direct  probe  of 
dynamical  processes  in  condensed  phases  [I-IH].  It  has  been  used  to  probe  relaxation  rates  and 
transport  phenomena  in  molecular  crystals  [2].  ionic  crystals  [3-6].  solutions  of  dye  molecules  [7.  S], 
conjugated  polymers  [9],  semiconductors  [1()-14|,  surfaces  [IS]  and  proteins  (16|.  Laser  gratings  are 
important  for  holography  and  the  elimination  of  beam  distortions  in  random  media  (phase  conjugation) 
[17.  18].  A  simple  phenomenological  description  of  the  TG  can  be  obtained  as  follows;  When  two  laser 
tields  with  wave  vectors  k.  and  k,  simultaneously  interact  with  the  material  svstem.  they  form  an 
interference  pattern  with  wave  vector  k,  -  A,.  Consequently  some  material  property  (e.g.  excitons. 
electron-hole  pairs,  space  charge,  temperature,  density,  etc.)  is  created  and  modulated  by  the  same 
wave  vector.  If  the  optical  dielectric  function  depends  on  this  property,  it  will  be  spatially  modulated: 
Af(r)  =  A  cos((A,  -  AO  •  r|.  When  a  third  beam  with  wave  vector  A,  is  now  scattered  by  the  sample,  it 
will  undergo  diffraction  resulting  in  a  signal  at  wave  vector  A^  =  A,  =  (A,  -  AO-  By  varying  the  delay 
between  the  initial  and  the  probe  beams  we  can  follow  the  decay  of  the  grating  amplitude  and. 
consequently,  the  underlying  motions  of  the  elementarv  excitations,  on  a  controlled  lengthscale 
determined  by  the  grating  wave  vector  A,  =k,  -  A,  (rtg.  I).  For  off-resonance  conditions  \e  is  purely 
real  (phase  grating)  whereas  for  resonant  excitation  it  has  an  imaginary  part  (population  grating). 
Phenomenological  treatments  iif  grating  experiments  have  been  developed  and  widely  used  for  the 
interpretation  of  numerous  experiments.  These  consist  of  identifying  the  nature  of  the  relevant 
dynamical  variables  produced  by  the  laser  grating  (excitation  density,  free  charge  carriers,  temperature, 
mass  density,  etc.)  and  using  macroscopic  transport  equations  to  describe  their  time  evolution.  .An 
excellent  review  of  this  level  of  description  is  found  in  Eichler  et  al.  [1]. 

signal 


k5=  ki  -  k2 -t-  k3 


Fig  I  Tvpical  transient  grating  setup  Fun  exeitation  beams  crossed  under  ancle  Oereate  a  crating  in  the  sample  with  wave  sector  k  -  k.  After  a 
sariable  delav  r,  the  grating  is  probed  b\  a  third  pulse  k  .  resultinc  in  a  nonlinear  i  ■  ditlracted  '  I  signal  M  k  -  k  k  .  -  k 


4 


/  Kn(ffMt'r  unii  S  Mukunwt.  I  riinsicni  ■^riutni!\.  !*i>tunion  :•  ’u^fi^incur  ot^tus 

Alternatively.  IG  speetroseopy  ean  he  viewecTas  a  particular  techniuue  ot  [unilinear  uptics  which 
belongs  to  a  broader  tamily  ot  processes  known  as  lonr-wavc  inixim’  |1.  1  our-waive  mixine  is  a 

technique  in  which  three  incoming  beams  with  wave  vectors  k..  k-.  k,  interact  with  a  nonlinear  medium 
and  generate  a  coherent  signal  at  a  wave  vector  which  is  anv  combination  of  the  wave  vectors 
-k^.-k..-k..  The  tour-wave  mixing  signal  is  usually  calculated  using  the  standard  and  well 
established  systematic  machinery  ot  nonlinear  optics.  In  this  picture  we  expand  the.  optical  polarization 
in  powers  ot  the  average  electric  held  ( the  Maxwell  Held ) 

P  =  \  '  E  +  \  '  - EE  -  \  EEE - 

Four-wave  mixing  is  related  to  the  nonlinear  susceptibilitv  \  The  tour-wave  mixinti  point  ot  view 
allows  us  to  put  the  TG  technique  in  a  broader  context,  and  to  explore  its  relationship  to  other 
spectroscopies  (photon  echo,  pump-probe. .CARS,  etc.)  |iy-27|.  For  example,  it  has  been  shown  that 
the  information  obtained  trom  the  FG  experiment  ean  in  principle  also  be  obtained  using  a  trequenev- 
domain  technique,  degenerate  tour-wave  mixing  (D4WM)  [28-35].  In  this  variant,  three  staiionarx 
incident  waves  with  wave  vectors  k,.  k.  and  k,  interact  in  the  svstem  and  the  scattered  wave  at 
-  1^:  ~  k,)  \s  detected.  It  has  been  shown  that  resonances  in  the  signal  as  a  function  ot  to,  -  lo, 
are.  under  some  very  general  conditions,  equal  to  the  Fourier  transform  ot  the  TG  signal  [3b].  The 
possibility  of  exeiton  localization  [37-30].  which  is  the  analogue  ot  the  .Anderson  electron  localization, 
could  also  be  probed  ideally  by  the  grating  technique  [4()[. 

In  this  review  we  develop  a  fully  microscopic  tnimework  tor  the  calculation  ot  four-wave  mixing 
processes  in  condensed  phases,  and  use  it  to  analyze  FG  spectroscopy  in  molecular  crvstals.  Molecular 
crystals  at  low  temperatures  seem  ideal  systems  (or  the  application  ot  the  grating  techniques.  In  the.se 
>ystems.  optical  excitation  creates  elementary  excitations.  Frenkel  excitons  [41 — tb]  ()r  charge  transfer 
excitons  [.■'()[.  which  are  well  understood.  Such  systems  seem  suthcientlv  simple  to  allow  a  rigorous 
microscopic  treatment  [5()-56|.  It  was  anticipated  that  the  grating  technique  would  directiv  probe  the 
exciton  motion.  Fhat  motion  is  expected  to  be  incoherent  at  high  temperatures  ;md  to  graduallv  become 
coherent  as  the  tempeiature  is  lowered.  Experiments  performed  on  anthracene  crvstals  [2b. c[  showed 
however  no  evidence  of  coherent  (nondiftusive)  exeiton  motion.  Instead  thev  showed  incoherent  motion 
with  a  very  large  diffusion  coetheient.  These  observations  were  interpreted  bv  .Agranovich  et  al.  [51.  52] 
in  terms  of  diffusion  ot  poUmions  [?7-b2].  which  are  quasiparticles  representing  the  correlated 
polarization  and  radiation  held  degrees  ot  freedom.  Other  evidence  related  to  polariton  dvnamics  in 
organic  crystals  was  obtained  by  .Small  and  coworkers  [b(l|  who  measured  second  harmonic  generation 
and  two-photon  Huorescence  in  naphthalene.  Fheir  measurements  strongi\  suggest  the  importance  ot 
polariton  (rather  than  exeiton)  scattering,  and  cannot  be  accounted  lor  using  standard  (nonretarded) 
expressions  tor  the  nonlinear  susceptibilities.  Polariton  ettects  have  also  been  measured  bv  a  varietv  of 
other  nonlinear  optical  techniques  [63-65 [. 

Fhe  incorporation  of  polariton  effects  in  the  theory  of  nonlinear  susceptibilities  is  not  straightforward 
(.''1.  .''2.  66-69|.  Traditionally  the  optical  susceptibilities  are  viewed  as  purelv  material  quantities,  and 
calculated  using  summations  over  eigenstates  of  the  molecular  (unretarded)  Hamiltonian.  Consequently 
V  depends  on  exeiton  resonances  and  dipole  matrix  elements  and  does  not  depend  on  retarded 
interactions  [21].  The  calculation  ot  any  nonlinear  optical  signal  is  then  convenientlv  divided  into  two 
steps;  we  hrst  calculate  susceptibilities  in  the  absence  of  the  radiation  held  and  next,  we  use  the 
susceptibilities  in  the  Maxwell  equations,  thereby  introducing  the  retarded  interactions  on  a  macro¬ 
scopic  level.  This  is  the  conventional  formulation  of  nonlinear  optics  developed  bv  Bloembergen 
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1 which  In'IJs  as  lone  as  ihc  ladialu'n  and  inaticr  dcercc"  o!  Irccdom  arc  wcaklv  coupled  le  e  . 
tor  i)tt-rcsonant  processes  or  when  exciton  dephasine  o  lareei  i  nder  these  eondiiions  the  eleineniarx 
excitations  are  excitons  and  photons,  and  it  o  isossihle  to  neelecl  correlations  between  the  radiation 
field  I'ld  the  material  polarization  and  to  detine  susceptibilities.  W  hen  the  radiation  and  matter  modes 
are  stronglv  coupled,  then  retarded  interactions  need  to  be  incorporated  more  microscopicalK  and  the 
optical  response  reflects  the  dvnamics  of  polaritons  it  is  then  ditticult  to  calculate  suscep¬ 

tibilities.  Instead,  the  suinid  is  expanded  directlv  in  terms  o|  the  cxicrtuil  field  rather  than  the  .\/iaut'// 
field  E.  It  should  be  stressed  that  the  linear  susceptibihtx  |hS.  .s(i|  \  land  the  dielectric  tunctum  r  I  can 
alwa\s  be  unambieuousK  defined,  whether  or  not  polariton  etteets  are  important.  ()ptieal  nonlmeanties 
require,  howexer.  a  different  approach  in  both  situations. 

One  i)t  the  primarv  eoals  ot  the  present  rexiew  !■'  to  dexcKip  a  systematic  tormulation  in  which 
polariton  and  exciton  dxnamics  can  be  accounted  ti'r  m  a  unified  framework.  In  addition,  we  discuss 
■-everal  approximate  schemes  tor  the  calculation  ot  nonlinear  optical  response.  Ihc  simplest  is  the 
anharmonic  oscillator  moilel  tor  the  polarization,  proposed  b\  Bloembcreen  |l‘i|.  Moreover,  manx- 
bodv  etteets  are  often  handled  using  a  mean-field  thex'rx  (the  U>cal-tield  approximation)  pO.  l|.  Czur 
general  formalism  reduces  to  these  common  procedures  xxhen  spccitic  approximations  are  made  and 
therefore  provides  useful  msieht  into  their  limitations. 

Ihe  outline  (d  this  rexiexx  is  as  tolloxxs  In  section  1  we  introduce  the  general  model  Hamiltonian  tor 
.1  molecular  crvstal  which  inckkfes  excitons.  piumons.  .md  photons  and  their  coupling  Ihe  multipolar 
{  ^  •  l)\  form  of  the  exciton-phonm  couplmu  is  adopted.  Ihc  entire  rexiew  is  dexixted  to  developing 
approximate  schemes  for  calculating  the  dvnamics  described  bx  this  Hamiltonian  and  the  Heisenberg 
equations  (eqs.  Z.d)  for  xarious  limiting  eases.  In  section  we  focus  on  the  material  evolution  alone 
I  excitons  and  ptumonsi  bx  factorizing  out  the  electn'iiiagnetic  field,  treating  it  effectively  as  a  classical  (■ 
number.  This  is  the  iismil  iramexvork  m  which  linear  and  nonlinear  optical  susceptinilities  are  defined 
and  calculated.  In  section  .v  1  we  present  the  operatx't  equations  and  their  solution  to  linear  (xrder  in  the 
electric  lield  \  ,  The  material  equations  form  a  hicrarchv  w herein  smele-pariiele  operators,  such  as 

the  polarizatK'fi.  are  cmiplcxl  lo  operators  mxoixing  sueeessixeix  more  particles  ](>  lowest  I'rder  in  the 
liierarchx  we  retain  onix  sinele-partide  operators  ,ind  describe  tour-waxe  mixim:  and  the  nonlinear 
^usceptlbilltx  \  ’  xviihin  the  loeal-field  .ipproximation  Kection  .'  2).  Hie  incorporation  id  transport 

requires  the  next  1  txxii-particlei  level  of  the  liierarchx.  wnich  is  adilressed  in  '•ection  T.V  Ihe  degenerate 
four-xvaxe  mixing  sional.  which  is  absent  in  the  smgie-particle  lexel.  noxx  ^hows  up.  and  dephasmg- 
mduced  resonances  are  obtained.  In  section  .'v4  we  discuss  the  transient  eratim:.  which  is  the 
■'ine-domam  .inah'eue  oi  oecenerale  'our-waxe  niixine.  !  In-'  o  onix  vtor.e  'xiiliin  I'le  'xxo-particle 
description,  as  no  transport  exists  on  the  smeie-partieic  lexel.  In  section  .v.''  we  ■'how  how  the 
wxo-partide  eouations  of  the  matter  can  be  related  to  the  Boltzmann  aiui  ibiC  diffusion  euuations.  We 
iurther  introduce  the  strone-eoilision  and  the  Haken-btrobl  mooels  tor  exciton-phonon  scattering.  I  his 
concludes  our  xliscussion  id  exciton  transport  and  dxnamics. 

In  section  4  we  turn  to  the  more  ueneral  and  complex  problem  xxhen  the  photon  xariables  are 
^tronelx  correlated  with  the  material  and  cannot  be  f.ietori/ed  In  this  ease,  we  lornuilate  nonlinear 
optics  xerx  diffcrentlx.  bx  usine  polaritons  and  axoidine  the  calculation  ot  snseeptibilitics.  In  section  4. 1 
we  introduce  the  polariton  transformation  which  exacilx  soixes  the  linear  optics.  In  section  4.2  we 
dexelop  equations  ot  motion  at  the  txvo-partide  lexel  i.malogous  to  section  .v.'l.  In  sections  4..''  and  !.  f 
xe  applx  these  equations  to  the  transient  grating  .md  its  frequency-domain  analogue  idegenerate 
tour-xvave  mixinc).  respectivelv.  detailed  discussion  of  xarious  limiting  cases  is  given.  Finallv.  m 
section  xvc  present  concluding  remarks  and  summarize  our  results. 


In  this  scctuni  uc  present  nur  model  system,  the  Hamiltonian,  and  the  hasie  equations  oi  inoiioii  iluii 
are  used  in  sections  and  4  to  calculate  optical  signals.  Wc  arc  interested  in  the  nonlinear  optical 
response  ot  molecular  crvstals  and  restrict  the  model  to  a  lattice  ot  polarizable  (nonpolar)  t\\o-ie\el 
molecules  with  transition  trequency  ii  (one  molecule  per  unit  cell).  The  total  number  ot  molecules  m 
the  crystal  is  denoted  bv  .\.  Tven  though  optical  experiments  are  usually  carried  out  on  thin  er\'-ial 
>labs.  \ye  will  assume  that  the  sample  is  still  thick  enough  ii>  invoke  translational  inv  ariance  in  all  hiiiiee 
directions,  hir  example,  we  will  assume  that  the  elementary  excitations  ot  the  crystal  lexcitoii''  or 
polaritons)  arc  well-approximated  by  those  tor  an  intinite  crystal.  Retinements  related  to  the  ire.itment 
ot  multilevel  molecules,  the  occurrence  ot  more  than  one  molecule  per  unit  cell,  or  the  explicit 
treatment  ot  tinite  geometries,  are  possible,  but  not  essential  teir  the  main  obiective  ot  this  paper  and 
would  obscure  the  phvsics  bv  adding  notational  and  algebraic  complexity. 

.-\n\  microscopic  thcorv  ot  (nonlinear)  optical  response  should  start  with  the  choice  ot  a  Hamiltonian 
sampling  the  material  system  to  the  radiation  field.  The  two  best-known  choices  are  the  multipolar 
{ fi-  D)  Hamiltonian  and  the  minimal  coupling  {p‘  A]  Htimiltonian.  which  are  related  bv  a  canonical 
trtmstormatum  1^2.  Here,  we  will  not  elaborate  on  the  controversial  issue  ot  the  exact  equivalence 
ot  both  choices;  an  extensive  literature  exists  which  deals  with  this  problem  (74|.  Suttice  it  to  suv  that 
discrepancies  between  the  results  of  both  pictures  are  ultimately  due  to  the  fact  that  approximations 
(which  are  unavoidable  m  a  practical  calculation)  affect  them  m  a  different  vvav.  We  will  use  the 
multipolar  Hamiltonian  tor  the  two  following  reasons;  (i)  This  choice  allows  for  a  straighttorward 
connection  tvi  the  literiiture  id  the  Bloch  equations,  which  are  used  to  describe  nonlinear  exciied-^tate 
vivnamics  v>t  isvilated  molecules  (22.  75|.  (ii)  An  easy  connection  to  the  popular  local-tield  apprv'ach  m 
condensed  systems  is  possible  within  this  Hamiltonian  176).  .A  drawback  iif  the  multipolar  Hamiltonian 
IS  that  it  does  not  explicitly  contain  intermolecular  interactions;  these  are  instead  carried  bv  exchange  ot 
photvms  between  the  molecules  [77],  The  interactions  can  be  recovered  bv  elimination  of  the  radiation 
held  [7S|  or  by  a  procedure  presented  previously  by  us.  that  keeps  the  radiation  held  as  a  vlcgrcc  ot 
freedom  (7h|.  The  latter  procedure  will  be  used  below. 

In  the  dipole  approximation,  the  explicit  torm  ot  the  multipolar  Hamiltonian  h/r  our  system  reiid'' 

r-’l: 

H  'In...  I  /’(D-/;  (r)dr-2-^  |  T,.,(r)! '  dr  *  T  i2.1i 

IthroughiHit  this  paper,  operators  are  indicated  by  a  caret  ((H;  the  ''amc  svmbvd  withvnit  card  vienotc'' 
!he  expectation  value  Oil)  d.>(/))|.  In  eq.  (2.1)  is  the  Hamiltonian  of  the  isolated  molecule  ni 
and  H  I  IS  the  contribution  from  the  tree  radiation  field.  In  secvmd  quantization,  both  can  be  expre''''Cd 
using  creation  and  annihilation  operators. 

// .  ■  hilBB  2.2  ) 

n  H,  ~  ^1  w  d* '  •  1  2..'  I 

i  I 

Here.  /T„  ( /f,„ )  denotes  the  creation  (annihilation)  operator  for  an  excitation  vm  molecule  ni.  I  hcsc 
v'perators  commute  tor  ditterent  molecules,  where''  '  for  anv  simile  molecule  they  obev 
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rtic  last  equality  deiinos  the  nioleeuiar  population  operator,  uhieh.  eombinmu  eqs.  (2.4a)  and  i2.'‘ei. 
can  also  be  written  as 

\V„=1B:,:B.:.  (I.cei 

Evidently,  this  operator  has  eigenvalue  zero  in  the  ground  state  and  2  in  the  excited  state.  For  harmonic 
oscillators  (bosons),  cq.  (2.4a)  does  not  hold  and  VV'_,  as  dclincd  m  cq.  (2.4b)  vanishes  idcnticaib. 
rheretore.  neglecting  this  operator  tor  two-levei  molecules  is  reterred  to  as  the  Bose  approximation 
[42.  58.  bl).  In  eq.  (2.."^).  the  operators  and  create  and  annihilate  a  photon  with  wave  \cctor  k 
and  transverse  polarization  A  (  A  =  1.2).  respectively.  Thev  obe\  Bose  commutation  relations. 

V  I  =<>• 

and  they  commute  with  all  material  operators.  Furthermore,  k  siA'i  and  c  is  the  velocity  ot  light  in 
vacuum. 

Ihe  third  term  in  eq.  (2.1 )  gives  the  interaction  between  the  radiation  held  and  the  molecules,  /’(r)  is 
the  polarization  lield  in  the  medium,  which  in  the  dipole  approximation  mav  be  written  as 

Bir) -  r,,.]  .  (2.oa) 

Here  fx,.  denotes  the  total  dipole  operator  ot'  molecule  m  (position  r„j.  which  can  be  written  as 

=  fxjB,,,- BJ .  '  (2.0b) 

M,.,  ^he  transition  dipole  matrix  element  ot  molecule  m.  I)  (r)  denotes  the  transverse  part  oi  tlie 
electric  displacement  held  ;it  position  r  and  is  related  to  the  transverse  jxirt  ot  the  Maxwell  electrie-neld 
operator  Eir)  bv 

b  ir)  =  E  (r)  ~^-P  (r)  .  ,2.") 

We  stress  that  in  the  multipolar  Hamiltonian,  the  displacement  lield  /)  (r)  (and  not  the  Maxwell  tield 
E  ]  is  the  conjugate  momentum  ot  the  vector  potential  .-llr).  so  that  m  second  quantization  1)  ir)  is 
totally  expressed  in  terms  ot  nulidtion  creation  and  annihilation  operators  onlv.  Explicitly,  we  have  I'.'j 

■A  ir)  [Uj,  expliA-r)  +  exp(-iAT)]t'^,  .  i2s.i) 

,  .  V  /  2rr/iAc  . 

I)  (/•)-i/^l  — - — )  expliA ' r)  -  dj,  exp(-iA •  r)|e.,  .  i2.sbi 

with  V  the  quantization  volume,  taken  equal  to  the  crystal  volume  in  all  calculations,  and  r-jj  A  ^  i.  2) 
the  transverse  unit  polarization  vectors  belonging  to  the  wave  vector  A.  We  work  in  the  C'oukmib  gauge, 
so  that  the  longitudinal  part  ot  the  vector  potential  vanishes.  Of  course,  the  longitudinal  part  ol  the 
displacement  tield  also  vanishes,  because  there  are  no  free  charges  in  our  system.  The  transverse 
electric  tield  in  second  quantization  is  now  defined  bv  eq.  (2.7).  in  combination  with  eqs.  (2.b)  and 
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(2.8b).  In  general,  the  longitudinal  part  of  the  eleetrie  lield.  E  .  does  not  vanish,  but  is  fully  determined 
by  the  polarization  field  through  the  relation  E  =  -AttP  . 

The  fourth  term  in  eq.  (2.1 )  is  a  self-energy,  in  which  P,„(r)  is  the  transverse  polarization  lield  caused 
by  molecule  m  onlv.  For  two-level  molecules  this  term  is  an  infinite  constant  which  does  not  contribute 
to  the  evolution,  so  that  it  is  often  omitted  completely  from  the  Hamiltonian.  Finally,  the  last  term  in 
eq.  (2.1).  t.  represents  the  kinetic  energy  related  to  the  nuclear  motion  (phonons). 

To  calculate  optical  response,  coupled  equations  of  motion  for  the  expectation  values  of  dynamical 
variables  (e.g..  the  polarization  field)  are  ideal.  Such  equations  can  be  obtained  from  the  expectation 
values  of  Heisenberg  equations  of  motion,  supplemented  with  some  factorization  approximation  to 
truncate  the  coupling  to  increasingly  more  complicated  variables  (sections  3  and  4).  .A  great  advantage 
of  an  equation-of-motion  approach  is  that,  unlike  the  density-matrix  approach  in  the  Schrddinger 
picture  [19-23].  it  does  not  involve  the  tedious  calculation  of  the  eigenstates  of  the  total  system. 
Instead,  we  focus  directly  on  the  relevant  dynamical  variables  which  carry  the  information  necessary  tor 
a  given  measurement.  The  complete  dynamical  information  as  given  by  the  eigenstates  is  usually  too 
complex  and  not  required  for  most  applications  involving  complex  systems  in  condensed  phases. 

We  have  shown  earlier  [76]  that  the  HeiL'‘'iberg  equations  of  motion  for  material  operators  within 
the  multipolar  Hamiltonian,  can  be  written  in  a  form  which  explicitly  contains  the  instantaneous 
Coulomb  interactions  between  the  molecules,  and  the  interactions  between  the  molecular  dipoles  and 
the  transverse  Maxwell  field  E  (instead  of  O  ).  It  is  essential  that  in  the  derivation  of  these  equations, 
we  split  the  (transverse)  displacement  field  according  to  eq.  (2.7).  i.e..  in  terms  of  the  transverse  electric 
and  polarization  fields.  An  alternative  approach  that  uses  D  ^  D  =  E  +  AttP.  leads  to  equations  of  a 
very  different  form,  in  which  the  instantaneous  interactions  are  not  readily  recognized.  It  is  straightfor¬ 
ward  to  extend  the  derivation  of  ref.  [76]  to  the  present  situation  where  also  nuclear  motion  is  possible. 
For  an  arbitrary  operator  Q  (material,  radiation,  or  mixed),  we  find  (all  operators  taken  at  time  t) 

(\  \)dQ'dt=  LQ  . 

=  [^n,.„  ^  ^1-  -  3  2  (-21  (Tm)  *  (Tj  •  (21! 

^  rn 

-  i  2  [M,„-[0‘(r,J.  ()[.  * \D'ir„].  C'j  •  /i,J  • 


(2.9a) 


(2.9b) 


Here.  denotes  a  material  Hamiltonian  which  consists  of  three  well-known  parts. 


is  the  usual  Frenkel  exciton  Hamiltonian  [42[ 


^2'Ar,„J(fl,*5,J(  ft, -«,„). 


12.10) 


12.11) 


where  the  second  term  accounts  for  the  instantaneous  dipole-dipole  interactions  between  the  molecules 
in  their  equilibrium  positions  and  orientations  (the  prime  excludes  terms  with  ni  =  n  from  the 
summation).  We  have  defined  s  r  -  r„.  where,  from  now  on.  r,„  denotes  the  equilibrium  position 
of  molecule  m.  The  explicit  form  of  the  interaction  reads 


fiJ(r)  =  ft  -  (I  ir'  -  }rrlr')‘  fi  . 


(2.12) 


J.  Knoester  and  S.  Mukamel.  Transient  gratings,  four-wave  mixing  and  polariton  effects  in  nonlinear  optics 

with  fi  the  molecular  transition  dipole  in  the  equilibrium  configuration  (equal  for  all  molecules  on  the 
lattice).  Of  course,  short-range  interactions  (e.g.,  exchange),  which  have  been  omitted  from  the 
Hamiltonian  from  the  very  start  (eq.  2.1),  can  be  added  heuristically  to  Jir).  The  second  term  in  eq. 
(2.10)  is  the  energy  of  the  phonons,  which  are  treated  within  the  harmonic  approximation  [42,  79), 

.  (2.13) 

q.s 

Here,  (fe^J  is  the  creation  (annihilation)  operator  for  a  phonon  with  wave  vector  q  in  branch  s  and 
is  its  energy,  and  obey  Bose  commutation  relations  [cf.  eq.  (2.5)]. 

Finally,  the  last  contribution  to  is  the  exciton-phonon  interaction,  which  arises  from  the 
dependence  of  the  intermolecular  interactions  and  the  Van  der  Waals  shifts  on  the  displacements  of  the 
nuclear  coordinates  from  their  equilibrium  values.  To  lowest  (linear)  order  in  these  displacements,  we 
can  write  [42.  60a] 


1 

VN 


2  [F,(*.  q)  +  +  ft-,,)  - 


(2.14) 


where  and  are  the  exciton  annihilation  and  creation  operators  in  momentum  representation, 
respectively. 


s; 


1 

V'V 


Z5;„exp(iA:T„,). 


(2.15a.  b) 


Of  course,  these  operators  are  periodic  on  the  reciprocal  lattice  and  the  inverse  transformations  read 
2  fi*exp(iA:T,J  .  A, 2  exp(-iA: •  r,J  .  (2.15c.d) 

IBZ  IBZ 

where  the  k  summations  extend  over  the  first  Brillouin  zone  only.  F^[k,  q)  and  xXr)  complex 
coupling  constants  that  can  be  expressed  in  the  first  derivatives  of  the  intermolecular  interactions  and 
the  Van  der  Waals  shifts,  respectively,  with  respect  to  the  nuclear  displacements  [42].  They  obey  the 
symmetry  relations  F^k.  q)  =  Fl{k  +  q.  -q)  and  \Sq)  ~  \*i~q)  (the  asterisk  denotes  the  complex 
conjugate),  which  guarantees  that  //.^  is  Hermitian.  It  is  noteworthy  that  in  the  delocalized  representa¬ 
tion  of  eq.  (2.14).  F^(k,  q)  and  x^{q)  multiply  the  same  operators,  so  that  their  sum  may  be  replaced  by 
a  single  total  coupling  constant,  which  we  will  write  F,{k.  q).  In  the  standard  reference  [42]  for  this 
is  not  the  case,  as  a  result  of  improper  use  of  translational  symmetry.  This  completes  the  discussion  of 
in  eq.  (2.9b). 

All  other  contributions  to  this  equation  speak  for  themselves.  In  the  last  two  r.h.s.  terms,  the 
equilibrium  positions  and  orientations  for  the  molecules  are  implied,  in  agreement  with  the  common 
neglect  of  direct  photon-phonon  interactions  [60.  80].  Clearly,  for  a  purely  material  operator  <2,  eq. 
(2.9)  reduces  to  eq.  (12)  of  ref.  [76],  with  additional  contributions  due  to  exciton-phonon  interactions. 
Equation  (2.9)  is  the  basis  for  all  equations  of  motion  that  are  used  in  this  paper.  In  general,  eq.  (2.9) 
will  result  in  an  infinite  hierarchy  of  coupled  dynamical  equations  whereby  single-body  operators  are 
successively  coupled  to  more  complex  quantities.  Fortunately,  the  optical  response  to  electromagnetic 
fields  that  are  not  too  strong,  requires  the  explicit  introduction  of  only  few-particle  states.  This  allows  us 
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to  truncate  the  hierarchy  at  a  very  early  stage.  We  shall  demonstrate  in  the  coming  sections  how  a 
truncation  at  the  two-particle  or  even  the  single-particle  level  is  adequate  for  the  calculation  of  a  variety 
of  optical  measurements.  This  situation  is  formally  very  similar  to  the  zero-temperature  many-body 
theory  where  a  few  quasiparticles  dominate  the  dynamical  behavior  [44.  8I|. 

We  conclude  this  section  by  giving  the  definition  for  the  spatial  and  temporal  Fourier  transforms  of 
an  arbitrary  function /(r.  r), 

f{k.(v)  =  jdr  J  d/ /(r, /)  exp(-iA  •  r  +  ioif) .  (2.16a) 

V 

The  inverse  transform  then  reads 


fir  t)  = 


J_  V 

IttV  ^ 


dct)  f{k.  <d)  exp(i^  •  r  -  icot) 


where  the  k  summation  extends  over  all  Brillouin  zones. 


(2. 16b) 


3.  Nonlinear  optical  response  of  excitons 

3.1.  Operator  equations  and  linear  optics 

In  this  section,  we  consider  the  exciton  theory  of  nonlinear  optical  response.  This  is  the  conventional 
approach,  in  which  the  electric  field  £*  is  treated  as  an  external  c-number.  Equivalently,  in  equations 
of  motion,  the  expectation  value  {QE~}.  with  Q  an  arbitrary  material  operator,  is  always  factored  into 
{Q){E^).  All  material  variables  can  then  be  expanded  in  powers  of  (£' ).  For  the  polarization  field, 
the  expansion  coefficients  define  the  susceptibilities  or  response  functions  (19-22).  In  combination  with 
the  Maxwell  equations,  the  susceptibilities  suffice  to  calculate  the  optical  signal.  In  this  approach,  the 
susceptibilities  are  completely  determined  by  the  evolution  of  the  isolated  material  system  with 
instantaneous  intermolecular  interactions,  i.e..  by  the  eigenstates  of  which  are  the  Coulomb 
excitons.  This  approach  does  not,  therefore,  account  for  the  fact  that  in  low-temperature  crystals 
strongly  mixed  coherent  combinations  of  photons  and  excitons  (polaritons  [57.  58] )  occur  as  eigen- 
modes;  polariton  effects  will  be  studied  in  section  4. 

The  first  step  in  the  exciton  theory  consists  of  deriving  equations  of  motion  for  the  exciton  operators 
and  Bl*  (with  k  in  the  first  Brillouin  zone)  from  eq.  (2.9).  Throughout  this  paper,  we  neglect 
Umklapp  contributions  to  £*.  i.e..  we  neglect  Fourier  components  E~{k.  t)  with  k  outside  the  first 
Brillouin  zone.  This  approximation  is  customary  (see,  e.g..  ref.  [57))  and  has  in  the  context  of  the 
multipolar  Hamiltonian  been  discussed  in  ref.  [76j.  Using  the  commutation  relation  (2.4b).  the 
following  compact  form  for  the  equations  of  motion  can  now  be  obtained  (all  operators  taken  at  time  t): 

[dlBA  i-n-j{k)-I(k)  -J{k)  \l  BA 

i  df  BiJ  ~  BiJ 

+  [2nph-'pL‘E^{k)  +  Mk)](  ) . 


(3.1) 
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II 


Here,  Jik)  is  the  lattice  Fourier  transform  of  the  intermolecular  interaction 

/(*)=  2  y(r,Jexp(-i*-r„,) .  (3.2) 

Because  the  lattice  model  described  in  section  2  is  centrosymmetric.  we  have  J{k)  =  J(-k).  Further¬ 
more.  I{k)  (=  .^'(-A:))  is  the  complex  exciton  self-energy  which  accounts  for  the  effect  of  phonon 
scattering  on  the  electronic  evolution.  In  appendix  A  we  discuss  a  procedure  to  calculate  I{k)  and 
obtain,  to  second  order  in  the  exciton-phonon  interaction,  the  well-known  expressions  [42| 


m^Mk)-ir{k), 

Mk)=-  -^,Pl\F^k.  ^)|'( 

ll  i\ 


(fhs) 


+  Or 


(3.3a) 


(3.3b) 


nt)  =  ^  2  - «,  - 4,)+  (« +  i)t«(A.,  -  tK  - « „)l .  (3-'0 


Here.  0^  is  the  frequency  of  the  Coulomb  exciton  with  wave  vector  k  (eq.  3.6).  denotes  the 

thermal  (Bose-Einstein)  occupation  number  of  phonons  with  wave  vector  q  in  branch  s.  and  P  stands 
for  the  principal  part.  Fik)  is  the  damping  rate  of  the  exciton  induced  by  the  phonon  bath;  A(it)  is  the 
phonon  induced  frequency  shift,  which  will  henceforth  be  neglected.  Fik)  plays  a  crucial  role  in  the 
occurrence  of  polariton  effects  (section  4.1). 

Returning  to  eq.  (3.1).  p  =  ,V/C  denotes  the  average  molecular  density  in  the  crystal  and.  hnally. 
M{k)  =  n^[J(k')\/Jj(B,.  +  pf,~'fi-E-(k').  mit-it')K  .  (3.4) 

k 

Here,  the  k'  summation  (as  from  now  on  all  wave  vector  sums  in  this  paper)  extends  over  the  first 
Brillouin  zone  only,  and  W{k  -  k')  is  the  lattice  Fourier  transform  of  the  population  operator. 


^  S  VV,„  exp[-i(A;-*')T„,] 


-  V  2  k  ^k  -k  '  • 

*' 


Let  us  first  consider  eq.  (3.1)  without  the  inhomogeneous  source  terms  multiplying  the  sector  (  ). 

This  defines  an  eigenvalue  problem  whose  solutions  are  the  annihilation  operator  for  the  Coulomb 
exciton  at  wave  vector  k  and  the  creation  operator  for  the  Coulomb  exciton  at  -k.  in  terms  of  B^  and 
[42].  The  Coulomb  exciton  frequency  is  determined  by  the  secular  equation  of  the  problem  and 
easily  found  to  be  -  i/’(it).  with 


=  {n[n  +  2j{k)]}' - .  (3.6) 

For  |y(Ac)|  /2.  which  is  almost  by  definition  the  case  in  molecular  crystals  [42.  4(S|.  this  yields 

=  {}  +  J(k).  This  is  known  as  the  Heitler-London  approximation  and  is  obtained  directly  if  one  does 
not  use  second  quantization.  In  this  approximation,  the  Coulomb  excitons  are  simply  created  (annihi¬ 
lated)  by  B^  (B^)  [42].  The  excitons  respond  to  the  source  term  in  eq.  (3.1).  which  contains  a  linear 
(~E")  and  a  nonlinear  (—.H)  contribution.  The  latter  is  the  source  of  nonlinear  optical  response  and 
vanishes  identically  for  harmonic  oscillators  [VT(A:)sOj.  We  note  that  in  the  present  paper  W{k)  is  the 
only  source  of  nonlinearities.  In  systems  with  multilevel  (and  polar  two-level)  molecules,  other  sources 
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arise  from  intermolecular  interaction  terms  that  are  cubic  and  quartic  in  the  molecular  exciton  creation 
and  annihilation  operators  [47.  59).  Such  terms  even  give  rise  to  nonlinearities  if  a  Bose  approximation 
is  applied,  which  is  usually  done  in  the  literature  considering  this  kind  of  nonlinearities.  In  reality,  both 
types  of  nonlinearities  will  occur. 

It  proves  useful  to  rewrite  eq.  (3.1)  by  introducing  the  variables 

P{k)^VN{B,  +  Bi,),  Vi:ik)sVyV(Bi^5lj.  (3.7a.  b) 

P{k)  equals,  up  to  a  factor  /i.  the  material  polarization  field.  Adding  and  subtracting  the  temporal 
Fourier  transforms  of  the  two  equations  contained  in  eq.  (3.1)  and  eliminating  V{k.  oj),  we  obtain 

{-[CJ  +  ir(it)]-  +  nl}Pik  w)  =  inph  'fi  -E  ik,  OJ)  +  .«(*,  «>) ,  (3.8a) 

Mik,  (o)  =  ^  [  diu' 'E  lJ(k')P(k\  a,')  -  •£*(*'.  w'),  W(k  -k'.oi-  Oi')],  .  (3.8b) 

ZiT  J 

Alternatively,  we  can  write 

„«(*,  CJ)  =  ,  fdoj'  I  doj"  EEin-aj"~  ir(r)][/}  +  co  -  co' -  co"  +  ir(k  -  k' -  k")\ 

{IttN)-  j  j  k'  r 

X  [J{k')P{k\  0)')  -{plh)ii-E-{k'.  a>'),  P{k".  io")P{k  -  k'  -  k'\  to  -  u)' -  o)")].  , 

(3.8c) 


where  W{k,  cj)  has  been  reexpressed  using  eq.  (3.5)  and  the  exact  relations 

B,{u))  =  [/]  +  01  +  ir{k)]P{k.  (o)I{2QV~N)  ,  (3.9a) 

Bl*(ai)  =  [n-o}-  ir{k)]P{k,  oi)/(2/2VA)  .  (3.9b) 

Here.  B*(ai)  denotes  the  temporal  Fourier  transform  of  Bl{t).  and  not  the  Hermitian  conjugate 
[6;t(oi)]'  of  B^(oj).  Of  course,  we  have  B'^{(o}  =  [Bj,(-a>)j'.  Relations  (3.9)  follow  easily  from  eqs.  (3.1) 
and  (3.7).  Transforming  eq.  (3.8a)  back  to  the  time  domain,  it  becomes  clear  that  we  have  replaced  two 
first-order  (in  time)  differential  equations  by  a  single  second-order  equation. 

To  conclude  this  subsection,  we  discuss  the  simple  case  of  linear  optics,  governed  by  the  first-order 
susceptibility.  We  first  note  that,  in  practice,  susceptibilities  are  defined  using  discrete  Fourier 
decompositions  for  the  expectation  values  of  the  fields,  instead  of  the  Fourier  transforms  eq.  (2.16)  [19. 
22.  68.  82).  The  electric  field  is  then  written  in  the  form 

E  ~(r.  t)  =  E  [Pj  e\p{'\k^  ■  r  -  m-t)  +  c.c.j .  (3.10) 

; 

where  /  labels  a  few  modes  which  are  essential  in  the  experiment,  oi,  ( >0)  and  are  related  by  the 
dispersion  relation  of  the  crystal.  The  Fourier  transform  eq.  (2.16)  for  this  field  reads 

E  {k.  oj)  =  IttV  E  [E;8,  ,  S{(o  -(0.)  +  E;*8,  _,  8((o  +  oi,)] . 


(3.11) 
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Similar  decompositions  are  possible  for,  e.g.,  the  polarization  field  P.  The  coefficients  appearing  in 
the  expansion  of  the  amplitudes  in  terms  of  powers  of  the  electric-field  amplitudes  E  * ,  are  now  the 
susceptibilities  [19,  22,  68,  82].  The  first-order  susceptibility  is  easily  obtained  from  eq.  (3.8a).  Namely, 
the  part  of  the  polarization  that  is  linear  in  the  electric  field  is  found  by  neglecting  the  nonlinearity 
Mik).  After  taking  expectation  values  and  substituting  eq.  (3.11)  and  its  analogue  for  the  polarization, 
we  obtain 


(3.12a) 


with  the  linear  susceptibility  tensor 


•d)/,  ,_£(*,  w)-l 

r  (k.(o)= - ;; - 

477 


2Qph  V/x 
-[(o  +  \r{k)Y- +  n;  ■ 


(3.12b) 


Here,  and  throughout  the  paper,  superscripts  in  parentheses  indicate  the  order  in  the  electric  field 
amplitudes.  e(fc,  w)  is  the  frequency  and  wave-vector  dependent  dielectric  tensor.  Equation  (3.12b)  is  a 
standard  result  (see,  e.g.,  ref.  [42]);  for  crystals  containing  molecules  with  more  than  two  well- 
separated  levels,  eq.  (3.12)  gives  the  contribution  of  each  level  to  the  linear  susceptibility.  In  the 
remainder  of  this  section,  we  will  address  the  nonlinear  optical  response  of  the  crystal  by  including  the 
effect  of  Mik)  in  various  approximate  ways. 


3.2.  The  single-particle  description:  frequency-domain  four-wave  mixing 

The  first  nonlinear  optical  technique  that  we  discuss  is  frequency-domain  four-wave  mixing.  We 
consider  ?  situation  with  three  fundamental  fields  [/=  1.2.3  in  eq.  (3.10)]  and  are  interested  in  the 
signal  at  i,  w,)  =  (k,  -  k,  +  k,,  w,  -  <u,  -I-  wj.  To  lowest  order  in  the  field  amplitudes  Ef,  this  is 
determined  by  the  third-order  susceptibility,  which  is  defined  through 

P';'*  =  =  -  w,;  k,w,,  -k.  -  a;,.  k,u)^)\E;E:*E:  .  (3.13) 

where  is  the  discrete  Fourier  coefficient  of  the  polarization  field  with  wave  vector  k^  and  frequency 
to  third  order  in  the  electric  field  amplitudes.  We  note  that  is  the  lowest  nonlinearity  allowed  by 
the  present  model,  since  y'"'  vanishes  for  a  centrosymmetric  medium  [21]. 

To  evaluate  we  take  the  expectation  value  of  eq.  (3.8a)  with  (k.  tv)  =  (k..  wj.  The  nonlinear 
source  term  for  the  third-order  polarization  is  {./((A:,,  wj)'’’,  which  involves  the  expectation  values  of 
products  of  two  and  three  "single-particle"  operators  [{£‘)(PP)  and  (PPP)  in  eq.  (3.8c)].  For  these 
products,  new  equations  of  motion  must  be  derived,  which  will  involve  yet  higher-order  products  of 
operators,  etc.  In  order  to  truncate  the  thus  generated  hierarchy  of  equations  of  motion  [61],  a 
factorization  approximation  must  be  invoked  that  breaks  apart  expectation  values  of  products  of 
operators  into  products  of  expectation  values.  It  is  natural  to  start  with  the  simplest  possible  truncation, 
which  consists  of  factoring  {.it)  completely  into  single-particle  expectation  values.  A  typical  contribu¬ 
tion  to  M^^\k^.  oij  then  reads 

a(/(*,)p;"  -  pft- V  •  £()F1"*PV'  • 

The  first-order  polarizations  in  this  expression  are  obtained  from  the  linear  approximation  to  eq.  (3.8a) 
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and  straightforward  algebra  leads  to 


-k^  - 


~k^  -  (t)^,  ky(o^)  =  Aflp 


fifififi  [fi  +  012-  ir(fe;)][/2  +  oj,  +  ir(fci)] 
d{k^.(oJ  A{k^,(i)^)  Ai-k2,~(02) 


( ~  ^  permutations  of  y  =  1, 2, 3  , 


(3.14a) 


A{k,  (o)  =  -[u)  +  irik)f  +  nl. 


(3.14b) 


The  permutations  of  (/c,,  o), ),  (-fc^,  -w,),  and  (fcj,  0)3)  account  for  the  different  time  orderings  with 
which  the  electric  fields  can  interact  with  the  sample  [23].  Expression  (3.14a)  can  be  simplified  if  the 
three  following  conditions  hold:  (i)  |7(/c)|  O  for  all  k,  which  is  generally  the  case;  (ii)  |a>y  - 12^  ( > 

|i(^y)|  (off- resonance  fundamental  waves);'(iii)  |  so  that  the  rotating  wave  approximation 

(RWA)  [75]  can  be  applied  to  the  signal  wave.  We  then  obtain 


3 )  / 

X  (~ftj  -  <Uj; , -ik,  -  0^3,  ifcjaij) 

=  2p{ti^fi,i/ft^)[aj^  -  {},  +ir(k^)]-'[(o^  -  Q,^  +  ir{k,)]-'[co2  -  n.^  -irik.)]'' 

+  permutations  of  /  =  1, 2, 3  .  (3,15) 


We  note  that  contributions  which  are  anti-rotating  with  respect  to  the  fundamental  frequencies  are  still 
contained  in  this  expression.  They  are  hidden  in  the  permutations;  for  example,  interchanging  (Ik,,  &>,) 
and  {-k2,  -oj^)  results  in  two  anti-rotating  denominators. 

It  is  instructive  at  this  point  to  make  a  connection  to  the  (damped)  anharmonic  oscillator  picture, 
which  is  a  popular  way  to  think  about  nonlinear  response  [19,  22].  If  we  take  the  expectation  value  of 
eq.  (3.8)  in  the  single-particle  factorization  and  transform  back  to  the  time  domain,  we  find  [neglecting 
r{k)~  with  respect  to  f?*] 


P{k- 1)  +  2r{k)Pik.  t)  +  nlP{k,  t)  =  2nph~'fi’E^{k,  t)  +  J{(k,  t) .  (3.16a) 

0  =  -  ^  2  2  [p^‘ V  •  E^ik\  t)  -  J(k’)P(k\  r)]{[/2  -  ir(/k")]R(/k".  t)  -  iP(k'\  t)} 

X  {[/2  +  iPik  -k'  -  r)JP(lk  -k’  -  k",  t)  iP(jk  -k'-  k'\  t)}  .  (3.16b) 


Equation  (3.16a)  represents  a  set  of  damped  harmonic  oscillators,  coupled  by  anharmonic  driving 
terms.  The  oscillator  picture  can  be  pushed  even  further,  if  we  realize  that  all  intermolecular 
interactions  can  be  eliminated  from  eqs.  (3.16)  by  introducing  the  local  electric  field  through 


ti-E^ik,t)^pi‘E^{k,t)-ih lp)J{k)P{k,  t) .  (3. 17a) 

For  dipolar  interactions  in  the  continuum  (A— *0)  limit,  this  coincides  with  the  Lorentz  local  field  [19. 
22.  83] 


0  =  E{k,  t)  +  5  TrP{k,  t) . 


(3.17b) 
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Here  we  used  J{k)  =  Airph  'fi  •  {kklk'  -  1/3)-  /i  fonhe  dipole  sum  [84-86]  and  E'^  =  as  there 

are  no  free  charges  in  the  system  (divZ)  =  0).  The  equation  of  motion  (3.16a)  can  now  be  rewritten 

P(ifc,  t)  +  2r{k)P{k,  t)  +  n^Pik,  t)  =  2nph~-'fi  •  E^^ik,  t)  +  M{k.  t) .  (3. 18) 

In  r-space,  this  equation  takes  the  compact,  form 

P{r,  t)  +  2rP{r,  t)  +  O'Pir  t)  =  0  +  Mr-  0  -  (3- 19a) 


with  r  an  arbitrary  lattice  site  and 

Mir,  t)  =  -ipfilhn)  •  E^ir,  r)[(/3  -  ir)P(r,  iP(r,  /)][(/}  +  ir)P(r,  t)  +  iP(r,  f)] .  (3.19b) 

Here,  the  k  dependence  of  the  damping  is  neglected;  otherwise,  a  local  picture  is,  of  course,  impossible. 
The  polarization  at  every  site  behaves  like  an  oscillator  driven  by  the  local  electric  field  and  by 
anharmonic  (nonlinear)  “forces”.  It  is  noteworthy  that  the  anharmonicity  is  a  function  of  both  the 
“position”  (P)  and  the  “velocity”  (P)  of  the  oscillator;  in  heuristic  anharmonic  oscillator  models,  one 
usually  ^’5:«:umes  an  anharmonicity  in  the  position  only  [19,  22]. 

The  fact  that  the  single-particle  factorization  used  here  leads  to  a  local-field  description,  is  not 
surprising  and  agrees  with  the  more  general  conclusion  that  any  theory  that  uses  a  factorization  of 
single-molecule  variables  is  equivalent  to  the  local-field  approximation  [76].  Consequently,  it  must  be 
possible  to  write  the  susceptibility  eq.  (3.14a)  as  the  third-order  molecular  polarizability 
y(-aj^;(u,,  -(u^.ojj)  multiplied  by  appropriate  local-field  correction  factors  [19,  22].  We  check  this 
explicitly,  y  is  easily  obtained  from  eq.  (3.14a)  by  setting  7(Jt)  =  0  and  p  =  1.  We  then  find 


-  (t)/,k^(OJ,  -k.  -  (1)2,  kjojj) 

i(a>J  i(aj,)  ^("3)  ^ 

=  o  -  -  -  -  y\^ (x)  \  o)y .  “ *  tu-j ) 

^  di-k2,-(02)  '  '  ■ 


j(oj)  s  -(a;  +  ify  +  {}-  . 


(3.:0a) 

(3.:0b) 


The  first  three  numerators  in  eq.  (3.20a)  simply  cancel  the  molecular  resonances  in  y.  We  now  restrict 
our  treatment  for  simplicity  to  the  case  that  all  wave  vectors  are  perpendicular  to  the  dipoles.  Then,  for 
optical  wave  vectors  and  dipolar  interactions  [84-86], 


Jik)  =  -{4Trl3fi)pn- . 


1) 


Combining  this  with  eq.  (3.12b).  it  is  easily  found  that 

i(a>) /Aik,  (o)  =  ][£(*,  cj)  +  2],  (3.22) 

so  that  eq.  (3.20a)  is  indeed  of  the  familiar  local-field  form  [19,  22].  We  finally  note  that  cascading 
contributions  to  which  are  usually  found  in  a  local-field  approach  [21,  82,  87],  are  absent  here, 
because  the  second-order  polarizability  vanishes  for  nonpolar  two-level  molecules. 
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3.3.  The  two-particle  description  :  degenerate  four^ave  mixing 

In  this  section,  we  extend  our  study  of  frequency-domain  four-wave  mixing  by  relaxing  the 
single-particle  factorization.  Instead  of  factorizing  {M(k,  w))  in  the  expectation  value  of  eq.  (3.8a)  [or 
eq.  (3.1)]  completely,  we  only  factor  the  population  W{k-  k')  from  it.  We  thus  do  not  break  up  the 
population  into  single-particle  variables,  as.  we  did  in  section  3.2.  We  then  obtainirom  eq.  (3.1)  [with 

(1/i)  d5,(0/dr  =  -[/2  +  J{k)  -  ir(*)|B*(/)  -  mBl,{t)  +  ^  •£-(it.  i) 

+  ^  2  [J{k')P{k\ /)  -  pr  V  •£-(*’,  t)]W{k  -  k\  t) ,  (3.23) 

and  a  similar  equation  for  fi!*(f)  =  (£!*(/)). 

W(k  -  k\  t)  now  appears  as  a  new  variable.  However,  instead  of  pursuing  the  hierarchy  by  deriving 
an  equation  for  the  population  itself,  it  proves  more  useful  to  introduce  the  two-particle  variables 

Q(k, 

=  ^  2  (fl,^(t)B„(0)exp[-i*-(r,„  +  r„)/2-t-ip-(r,„-r„)].  (3.24a) 

m,n 

which  are  the  diagonal  elements  (/c  =  0)  and  coherences  (ikT^O)  of  the  exciton  density  matrix  in  the 
momentum  representation.  From  the  last  form  of  eq.  (3.24a).  it  is  clear  that  k  is  conjugate  to  the 
exciton  center  of  mass;  p  is  related  to  the  classical  exciton  momentum.  The  significance  of  the  k  and  p 
variables  may  be  clarified  by  switching  to  the  Wigner  representation  for  the  exciton  density  matrix.  This 
is  done  in  section  3.5.  where  we  also  make  the  connection  to  common  transport  equations  such  as  the 
Boltzmann  equation.  The  exciton  popuhtion  can  now  be  written  [cf.  eq.  (3.5)]. 

m.  0=^2  Q{k-  P- 1) .  (3.24b) 

As  next  step  in  the  hierarchy  we  now  consider  the  equation  of  motion  for  Q{k,  p,  t).  We  first 
concentrate  on  the  electronic  (coherent)  part  of  this  equation,  i.e..  without  accounting  for  phonons. 
This  is  obtained  from  (it,  s  p  -  k/2:  k.  =  p  +  kl2) 

(d[5;,(r)£J0]/df)  =  ([d£;(/)/dr]flJ0)  +  {B,^{t)dB,it)idt)  . 

and  eq.  (3.1)  without  the  self-energy  terms.  In  the  final  result  the  following  approximations  are  made: 
(i)  (£')  is  factored  from  all  other  variables  (exciton  theory!);  (ii)  we  neglect  terms  coupling  to 
variables  of  the  form  {Bl^{t)B'_^^{t))  and  (r)Bj  (r)),  which  is  equivalent  to  invoking  the 
Heitler-London  approximation  on  this  level  of  the  hierarchy;  (iii)  we  neglect  all  variables  which  are 
higher  than  bilinear  in  the  exciton  variables  (such  as  (VVB)).  as  they  eventually  result  in  contributions 
to  the  polarization  that  are  of  order  four  and  higher  in  the  electric  field  amplitudes.  The  leading  order 
for  four-wave  mixing  processes  is  three.  The  phonon  (incoherent)  contribution  to  Q{k,p.t)  is.  to 
second  order  in  the  exciton-phonon  interaction,  derived  in  appendix  A.  The  total  equation  of  motion 
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then  reads. 

T  ^  Qik,  p,t)^[J{p-  kll)-  Jip  +  kl2)\Q(k.  p.  0-E  -(*:  P-  p')Qik.  p'.t) 

I  Qi  p 

-(VN/fiV)[B^^,  ,{t)pi‘E-{-p  +  kiZ.t)-  B;,,  ,(i)ti-E  (p^ki2.i)].  (3.25) 

The  first  r.h.s.  term  describes  the  free-exciton  motion  and  the  last  term  represents  a  source  for 
two-particle  coherences  created  by  electric  fields  in  a  sample  in  which  a  polarization  already  exists.  The 
second  r.h.s.  term  in  eq.  (3.25)  is  due  to  the  phonon  bath,  where  the  complex  self-energy  matrix  has 
the  form 


lik:  p.  p')  =  [I{p  +  k/2)-  I*(p-  ki2)]8^^ 

+  I{p  +  kl2,  p  -  k/2.  p'  -  p)  -  p  -  kl2.  p  +  kl2.  p'  -  p) .  (3.26) 

Here,  I{k)  has  been  defined  in  eq.  (3.3)  and  I{k^.  k..  k^)  is  given  in  eq.  (A.  13).  Note  that  I(k\  p,  p  ) 
consists  of  the  independent  sum  of  the  self-energies  of  the  excitons  2  2  ^^at  make  up  the 

coherence  Q{k,  p,t)  (T, -processes)  [22,  75]  and  additional  contributions  (J)  due  to  the  correlated 
dynamics  of  these  excitons  {T*-  or  pure  dephasing  processes  [22.  75]).  Due  to  the  pure  dephasing 
contributions,  the  single-particle  factorization  in  section  3.2  breaks  down;  this  will  be  seen  explicitly  in 
the  result  eq.  (3.33).  It  is  important  to  observe  that  the  phonon  bath  only  induces  coupling  (scattering) 
between  coherences  with  different  p  values.  The  variable  k  is  conserved  in  eq.  (3.25).  which  is  a 
consequence  of  the  system's  translational  symmetry  and  the  thermal  average  performed  over  the 
phonon  bath  in  appendix  A.  We  also  note  from  eq.  (3.26)  that  for  ^;  =  0  the  self-energy  is  purely 
imaginary.  The  physical  explanation  for  this  is  that  a  diagonal  density  matrix  element  has  no  frequency 
associated  with  it.  and  hence  no  frequency  shifts  either.  Equations  (3.23)-(3.25)  form  a  closed  set  and 
govern  our  two-particle  description  of  excited  state  dynamics. 

To  study  frequency-domain  response,  we  apply  temporal  Fourier  transforms  and  introduce  the 
variable  P  defined  in  eq.  (3.7)  to  obtain 

{-[(u  +  ir(jt)]'  +  ni}P{k.  uj)  =  2nph~'fi •  e  ik. w) 

+  ^  I  da)'^\J(k')P(k'.  uj')-ph  'n-E  Ik'.  uj‘)\  ^  Z  0{k-k.  p.  w  -  cu'i .  (3.27a) 

277  J  *  ‘  P 

and,  using  eqs.  (3.9), 

-(oQ(k,  p.  w)  =  [y(  p  -  k/2)  -  y(  p  +  k/2)]Q(k.  p.io)-^  ^{k:  p.  p')Q{k,  p'.  a») 

P 

-  j  dw'  ([/2  +  co'  +  if(  p  +  ki2)\P{  p  -t-  k/2.  u)')pL-  E  (-p  +  kiZ.  u)  -  co') 

-  [n  -0)'  -  in-p  +  it/2)]P(-p  +  kl2,  cj')n  E  (p  +  k/2.  w-oj')}  .  (3.27b) 

These  equations  have  to  be  iterated  in  order  to  obtain  the  third-order  susceptibility. 
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Instead  of  treating  the  general  case,  however,  we  will  direct  our  attention  to  the  degenerate 
four-wave  mixing  (D4WM)  experiment  [28-35).  In  this  experiment,  two  fundamental  fields  exist. 
(A:,,  oj, )  and  (A:,,  wj  and  the  nonlinear  signal  at  =2a;,  -  cu.  is  observed  while  the 

frequency  difference  w,  -  w,  is  tuned  through  zero.  We  will  assume  that  .  w,.  and  2a;,  -  co.  are  all  far 
from  the  material  resonance  and  we  solely  concentrate  on  contributions  to  the  third-order  polarization 
that  have  possible  resonances  at  w,  -  w,.  Furthermore,  we  restrict  our  study  to  the  Haken-Strobl 
model  for  the  phonon-induced  self-energies  [54.  88).  Within  this  simple  model,  the  basic  equations  are 

ld<B,.(l)>/dl^=-!(r-y)(B, .(/)>,  (3,:8a) 

|d(B;(()B.(,))/dr|,.  =  -ini  -  d.J  T  yKB  J/)B.(0)  ■  (3.:8b) 


where  [•  •  jp,,  denotes  the  phonon  contribution.  These  equations  imply  for  the  self-energies 

I(k)^-if(k)  =  -i(r+y)/2  .  I(k:  p,  p’)=  -i(r+y)S^,  +  if/iV  .  (3.28c. d) 

The  parameters  F  and  y  are  usually  taken  real,  so  that  only  imaginary  self-energy  contributions  are 
included.  F  is  the  pure  dephasing  rate  and  y  represents  the  population  relaxation  rate.  The  main 
motivation  to  use  this  model  is  that  it  allows  for  analytical  results,  while  still  preserving  the  essential 
physical  aspects  related  to  pure  dephasing.  Previous  studies  of  D4WM  (using  response  theory  [36])  and 
of  transient  gratings  [40.  55]  in  crystals  have  also  used  the  Haken-Strobl  model. 

Given  our  microscopic  expressions  for  the  self-energies  derived  in  appendix  A.  we  can  in  principle 
improve  the  theory.  Yet  alternative  (Markovian  and  non-Markovian)  exciton  self-energy  models  have 
been  proposed  within  the  theory  of  optical  absorption  in  molecular  crystals  [89-91].  However,  both  the 
evaluation  of  more  realistic  self-energies  and  the  subsequent  solution  of  the  equations  of  motion  will 
involve  extensive  numerical  calculations.  Other,  probably  more  practical  improvements  over  the 
Haken-Strobl  model  lie  in  the  connection  between  the  equation  of  motion  (3.25)  and  the  Boltzmann 
equation  (section  3.5).  This  connection  enables  us  to  utilize  the  well-developed  methods  for  solving  the 
Boltzmann  equation  within  transport  theory. 

We  now  calculate  the  third-order  susceptibility  for  the  D4WM  setup.  Define  k^  =  k^-k.. 
P^  =  {k^  +  Ac,)/2  (g  stands  for  "grating")  and  let  Q^{p)  denote  the  component  at  wave  vector  k^  and 
frequency  w,  -  to.  in  the  discrete  Fourier  decomposition  of  the  two-particle  coherences.  From  eq. 
(3.27a)  we  then  obtain  for  the  third-order  polarization  at  the  signal  wave  vector  and  frequency  and  with 
a  possible  resonance  at  w,  -  to.. 


4/2p 


Mk,. 


(1) 


[J(k,)P\"-pfi  ^h-e;\2q[-\p). 


(3.29) 


has  already  been  .solved  in  section  3.1  and  from  eq.  (3.27b)  it  follows  that  the  for  different 

p  values  obey  the  coupled  equations 


w 


,  -  )(?';’( p)  =  \J(p-  kJ2)  -  y( p  +  kJ2)  +  dF+y  )]0';’(  p) 


-  i  ^  S  C'/’( p')  +  (Tf^^  + 


(3.30) 
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Here,  and  are  "source  terms"  that  originate  from  the  last  r.h.s.  term  of  eq.  (3.27b). 

Using  eq.  (il2b)  for  the  first-order  polarizations,  we  hnd 


o-R  =  -ph^-[ 


■,1  fi  +  (x)^  +  i(r  +  y)i2  n  +  (D-.  —  r  +  y ) i 2 


-Oi, ) 


fiE.fiE: 


=  -  0^+i{r  +y)/2]“'  -[w, -i(/V  y)/2]'}/A  •£,  M  ‘  . 

,(  O  -  io^-  '\ir+  y)l2  O  -  (ti.  +  \{r+  y)  2  \  ^  ^  ^ 

- ^e-Ein-Ei 


(3.31a) 


^  ph  ‘{[a;,  -r  +  [{F  y)/2]  '  -  (<a,  -r  -  i(/'-  y),2]  ‘}/x  •  E ,  /x  ‘  E:*  . 


(3.31b) 


where  the  approximations  hold  if  jyfilc,  )|.  iy(A-,  )| ‘^ /i.  Obviously.  and  are  the  rotating 

and  anti-rotating  sources  for  the  two-particle  coherences,  respectively,  and  they  transform  into  each 
other  if  (A,,  w, )  and  (-A,,  -u)-,)  are  interchanged.  In  section  3.2  the  same  interchange  of  rotating  and 
anti-rotating  terms  under  permutation  of  the  fundamental  fields  was  observed  (eq.  3.15). 

Equation  (3.30)  is  solved  in  appendix  B  and  yields 


y  s  ^  ^ _ ^  ) 

p  “  -iw,,  -  i|y(A,)  -  y(A, )1  +  (/' -^  y )  -iw.,  -  ifyiA, )  -  y(A,)l  +  (/'+ y) 

(  I  “  V  ^  ~  ^  ~  ^  ^  •  (3.32) 

.V  p  • 


where  w,,  =  <a,  -  ta,.  Combining  eqs.  (3.29).  (3.12b).  (3.31).  and  (3.32).  the  third-order  susceptibility 
is  easily  found.  To  simplify  the  result,  we  apply  the  RWA  with  respect  to  the  fundamental  fields 
(o-^  =0)  as  well  as  the  signal  field.  Furthermore,  we  use  w,  -  ^}^  ^UlA. )  (w,  is  off-resonance),  so 
thaty(A|)P',‘*  in  eq.  (3.29)  is  negligible  relative  lo  ph  'fi  ’  E ^  .  Finally,  we  use  the  approximate  identity 

(Up  /y.  -  +  \(r  +  y) 

I _  ^  J 

(u,  ,  -t-  J(k.  )  -  y(A, )  +  i(r  -t-  y ) 

This  relation  becomes  an  exact  identity  within  the  Heiticr-London  approximation  (HLA)  /ij  = 
/2  +  y(A).  Here,  we  do  not  obtain  an  exact  identity,  because  the  HLA  has  not  been  invoked  in  eq. 
(3.29)  for  the  polarization,  whereas  it  is  used  in  the  equation  for  the  two-particle  variables.  We  finally 
obtain 


'(- Ap  -  (i»  ;  A,  (U| .  -  A,  -  lu,.  A,  (u, ) 

=  2p{flfJ^f^f^h')l^t)^-^}^+]{^^y)i2]  '[(u,  - /i,  *  i(/'- y)  2]  jiu,  - /ij,  -  i(/  ^  y )  2| 
h 

^(l-  ^^\-iuj,,-iJ(p-k^.2)^U{p-k^^2)^  r-y\']  .  (3..V) 

P 

We  immediately  observe  that  this  result  is  the  product  of  the  rotating  part  of  eq.  (3.15).  which  was 
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obtained  in  the  single-particle  factorization,  and  a  dephasing-induced  factor  which  is  unity  if  I  ~  •*  In 
other  words,  in  the  presence  of  pure  dephasing,  the  single-particle  description  breaks  down. 

Loring  and  Mukamel  (36 j  were  the  first  to  calculate  the  D4WM  susceptibility  for  a  molecular  crystal, 
within  the  Haken-Strobl  model  for  dephaslng  and  using  Liouville-space  Green-function  techniques. 
Specific  application  was  made  for  a  one-dimensional  system  with  nearest-neighbor  interactions,  lor 
which  the  sum  over  the  momenta  p  in^.  -(3.33)  can  be  evaluated.  Their  analysis- of  the  D4WM  signal 
for  various  limiting  cases  is  also  possible  our  general  expression  (3.33).  The  signal  intensity 
S{k^,  u)^)  is.  within  the  slowly  varying  amplitude  approximation,  proportional  to  U*  ~  *i^i- 

-Jfe.  -  w,,  )|'.  Any  frequency  dependences  deriving  from  the  first  three  denominators  in  eq.  (3.33) 

may  be  neglected,  as  we  assumed  that  w, ,  o),  a^  ca,  are  off-resonance.  The  only  frequency  dependence 
in  S{k^,  a>J  then  emerges  from  the  last  factor  in  eq.  (3.33). 

Three  limiting  cases  are  now  of  special  int^st  [36].  First,  we  consider  non-interacting  molecules 
(y(it)  =  0]  and  find 

S(k,,u)Jy^\  +  r[r  +  2y)lico]^  +  y^).  (3-34) 

This  signal  has  a  Lorentzian  resonance  at  w,,  =  0,  whose  width  is  the  inverse  of  the  excited-state 
lifetime.  The  resonance  vanishes  in  the  absence  of  dephasing.  These  dephasing-induced  resonances 
have  been  observed  by  Bloembergen  et  al.  [31]  in  the  gas  phase  and  have  been  denoted  PIER4 
(pressure-induced  extra  resonances  in  four-wave  mixing).  In  fig.  2  we  show  the  variation  of  such  a 


Fig  2  PIER4  signal  characteristics  m  Na  as  a  function  of  buffer  gas  pressure  |.^Ib|,  (a)  Ratio  of  peak  height  to  nonresonant  signal;  |h)  resonance 
full  width;  (cl  integrated  intensity  of  resonant  signal. 
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resonance  in  Na  induced  by  the  buffer  gas  He.  (The’dephasing  rate  /'is  in  this  case  proportional  to  the 
He  gas  pressure.)  In  molecular  crystals  similar  resonances  have  been  observed  by  Hochstrasser  et  al. 
[32]  and  have  been  denoted  DICE  (dephasing-induced  coherent  emission).  An  example  is  displayed  in 
fig.  3a  and  3b.  This  resonance  is  similar  to  the  one  discussed  here,  except  that  it  occurs  at  w,,  equal  to 
some  vibrational  frequency  rather  than  w,,  =  0.  This  is  a  Raman  transition  which  is  totally  analogous  to 
the  D4WM  resonance.  In  Hochstrasser's  experiment,  tu,,  =  -747  cm"'  represents  a- vibrational  mode 
and  in  Bloembergen’s  case  oi,,  =  17 cm"'  is  the  splitting  of  the  NaD  lines.  The  unique  and  surprising 
aspect  of  these  resonances  is  that  usually  dephasing  results  in  loss  of  coherence  and  line  broadening, 
whereas  here,  it  induces  new  sharp  resonances  as  w,,  is  varied.  The  reason  is  a  delicate  interference  of 
various  terms  contributing  to  which  exactly  cancel  in  the  absence  of  dephasing.  The  addition  of 
dephasing  eliminates  this  cancellation  and  results  in  the  new  resonance  [34.  23]. 

As  a  second  special  case,  we  consider  molecules  with  arbitrary  interactions  in  the  absence  of 
dephasing.  For  /"  =  0.  the  Haken-Strobl  model  describes  coherent  exciton  motion  on  the  lattice,  and 
from  eq.  (3.33)  it  is  clear  that  the  D4WM  signal  exhibits  no  resonance  as  a  function  of  a>,2  in  this  limit. 

We  finally  discuss  the  case  of  finite  interactions  in  the  strong  dephasing  (incoherent  or  diffusive) 
limit,  defined  by  r>\J{p-  kJ2)  -  J{p  +  kJ2)\  (for  all  p)  and  r>y.  In  this  limit,  the  Haken-Strobl 
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Fig.  3(a).  Stationary  four-wave  mixing  (coherent  Raman)  spectrum  of  a  doped  crystal  (pentacene  in  benzoic  acid)  as  a  function  of  temperature 
|32b|.  The  portion  of  the  spectrum  near  01,,  = '47  cm  '  is  shown  on  an  expanded  vertical  scale  Note  that  the  dephasing-induced  band  at 
iv.i  =  747 cm  '  grows  relative  to  the  peak  at  (j,,  =  755  cm  '  with  increasing  temperature.  ( b)  Growth  of  the  pure  dephasing  rate  /  ( T)  as  a  function 
of  temperature  (32b].  as  derived  from  experimental  measurements  of  the  intensity  of  the  dephasing-mduced  band  at  ai,,  =  747  cm  The  solid  line  is 
a  fit  to  an  Arrhenius  form. 
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model  describes  diffusive  exciton  motion.  In  appendix  C  it  is  shown  that  for  arbitrary  interactions  and 
dimensionality  we  now  have 


5(*,,  wjacl  +  — 


f(f  +  2y) 


co-;,  +  [y  +  {k,-k,lDj^' 
with  the  exciton  diffusion  constant  (tensor)  defiffed  by 


(3.35) 


(3.36) 


In  the  last  step  of  eq.  (3.36)  we  have  used  eq.  (3.2).  In  the  incoherent  limit,  the  D4WM  signal  shows  a 
Lorentzian  resonance  with  a  width  that  is  the  sum  of  the  inverse  excited-state  lifetime  (y)  and  a 
contribution  from  exciton  motion  (fJj^g)-  The  same  result  has  been  obtained  for  a  disordered  medium 
[82].  In  general,  is  a  tensor  depending  on  the  magnitude  and  the  direction  of  k^-,  if,  however,  the 
range  of  interactions  is  small  compared  to  the  characteristic  length  scale  \k^\~^  of  the  experiment,  eq. 
(3.36)  reduces  for  an  isotropic  cf-dimensional  system  to 

=  (3.37) 

dl  fn 


Alternatively,  eq.  (3.36)  may  be  written  as 

2  [*e  •  Wc(p)]' .  '  (3-38a) 

with  v^{p)  =  V^J{p),  the  exciton  group  velocity  at  wave  vector  p.  For  an  isotropic  </-dimensional 
system,  this  yields 


with  F.  an  effective  exciton  group  velocity. 


(3.38b) 


3.4.  Time-domain  four-wave  mixing:  transient  grating  and  exciton  dynamics 

We  now  turn  to  the  study  of  transient  grating  (TG)  experiments  within  the  exciton  theory.  The 
following  typical  setup  is  considered  (fig.  1).  At  time  f  =  0,  two  short  excitation  pulses,  (fc,,  w, )  and 
(fc,,  w,),  crossed  under  an  angle  0  interfere  in  the  sample  and  create  an  excitonic  grating.  The  decay  of 
the  grating  as  a  result  of  dephasing  and  population  relaxation  is  monitored  by  applying  a  probe  pulse. 
(*,,  wj,  at  f  =  T.  The  observable  is  the  time-integrated  intensity  of  the  nonlinear  (“diffracted”)  signal 
with  wave  vector  /c,  =  it,  -  it,  +  jfc,  and  frequency  =  w,  -  w,  -i-  w,  as  a  function  of  the  pump-probe 
delay  t.  The  electric  fields  now  take  the  form  of  eq.  (3.10),  except  that  £■"  (;  =  1, 2, 3)  is  replaced  by 
the  pulse  envelope  E^(t).  For  simplicity,  we  will  consider  square  pulses  with  polarization  parallel  to  the 
molecular  transition  dipoles,  with  amplitudes  E^,  and  with  duration  r-.  All  pulses  are  long  compared  to 
an  optical  cycle  {u)^r-  ^1),  but  short  (delta  pulses)  on  the  dynamical  time  scales  in  the  sample. 
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To  describe  this  experiment  microscopically,  we  use  our  equations  of  motion  in  the  two-particle 
factorization  [eqs.  (3.23)-(3.25)].  Working  with  the  polarization  as  variable,  eq.  (3.23)  takes  the  same 


form  as  eq.  (3.16a),  but  now  with 

0  =  2/2  S  lAk')P(k\  t)  -  ph-'n>E^ik\  r)]W(ifc  -  k',  t) .  (3.39) 

(fhird)  order  in  the  pulse  amplitudes,  the  observed  polarization  P{k^,t)  is  sourced  by 
■  U  (*,,  r)  in  eq.  (3.16a),  and  the  only  important  contribution  to  t)  in  eq.  (3.39)  arises  from 

the  term  k  =  k^.  (In  principle,  third-order  contributions  are  also  obtained  from  k'  =  ifc,  or  k'  =  ~k2,  but 
these  are  very  small  due  to  negligible  temporal  overlap  of  pump  and  probe  pulses.)  We  assume  an 
off-resonance  probe  pulse  ([wj  - /2|t3  ^  1),  so  that  the  polarization  envelope  associated  with  it 
instantaneously  follows  the  electric  field  [92], 

^3'’(0  =  ?”(*3>  ‘Oa)  •  E^{t) .  (3.40a) 

The  envelope  of  the  nonlinear  source  term  then  obeys 

(3,40b) 

Within  the  slowly  varying  envelope  approximation  [22],  the  observable  now  reads 

5(T)a|  -fc-r)!’.  (3.41a) 

r 

with  d,  the  probe  pulse  area  defined  by 

(3.41b) 


In  eq.  (3.41a)  we  used  the  fact  that  E;{i)  is  a  delta  function  on  the  dynamical  time  scale,  so  that 
W{k^  -  k^,  t)  does  not  vary  appreciably  during  the  pulse.  According  to  eq.  (3.41).  the  TG  is  determined 
by  the  Fourier  component  of  the  exciton  population  at  the  grating  wave  vector  k^  =  it,  -  ifcj  and  at  time 
^  result  agrees  with  the  usual  diffraction  picture  [1.  2a.  55.  93]  and  was  first  derived 

microscopically  by  Loring  and  Mukamel  [56],  using  response  theory.  These  authors  also  pointed  out 
limitations  of  the  diffraction  picture  by  showing  that  in  general  the  temporal  profile  of  the  signal  does 
not  follow  that  of  the  probe  pulse. 

The  evolution  of  the  population  can  be  found  from  eqs.  (3.24)  and  (3.25).  During  the  pump-probe 
delay,  the  last  term  in  eq.  (3.25).  containing  the  electric  field  amplitudes,  is,  of  course,  absent.  As  in 
the  previous  section,  we  restrict  ourselves  to  the  Haken-Strobl  model  for  the  phonon-induced 
self-energies.  If  we  define  the  Laplace  transform  of  a  time-dependent  function  /(/)  by 

X 

/(s)  =  J  df/(i)exp(-sf) .  (3.42) 

0 


then  eq.  (3.25)  translates  into 
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p.  s)  =  [U{ p  -  kJ2)  -U{p  +  kJ2)-  (r+y)]Q^-%.  p.  s) 


+  p\s)  +  p.t  =  0) . 


(3.43) 


This  equation  has  the  same  form  as  eq.  (3.30),  governing  the  D4WM  experiment,  except  that  -io)  is 
replaced  by  s  and  the  source  term  is  now  determined  by  the  grating's  initial  condition  p,  t  =  0) 

right  after  the  pump  pulses. 

Before  solving  eq.  (3.43),  we  first  determine  these  initial  conditions.  As  the  pump  pulses  are  short  on 
the  dynamical  time  scale,  the  phonon  bath  has  ho  time  to  establish  correlations  between  the  ket  and  bra 
sides  of  the  exciton  density  matrix  (the  crystal  State  is  still  a  pure  state),  so  that  immediately  following 
the  pump  pulses  we  may  use  the  factorization 


0'^’(A,,p,t  =  O)=(B;  _(0))"><fl_(0))‘ 


(3.44a) 


The  first-order  exciton  amplitudes  can  be  obtained  by  integrating  the  linear  part  of  eq.  (3.23)  during  the 
pump  pulses.  In  the  RWA,  assuming  |wy-42|T^<^l  {j=  1.2),  and  neglecting  intermolecular  interac¬ 
tions  during  the  pulses,  we  find  [56] 


(3.44b) 


and  (6^(0))^'’  -  [{6*(0))‘'*]*.  Here  the  t?,  (/=  1,2)  denote  the  pump  pulse  areas  defined  in  analogy 
with  eq.  (3.41b).  Thus 


Q''-\k^,p,t=^Q)==pVd,d*8  , 


(3.44c) 


with  Pg  —  (^,  -I-  k.)l2.  as  before.  As  mf'st  theories  for  the  TG  signal  work  within  the  site  representation, 
it  is  useful  to  mention  that,  using  eq.  (2.15),  the  present  initial  condition  (a  coherence  of  excitons  at 
wave  vectors  fc,  and  it,)  translates  into 


( =  1^1'  cos(A:^  •  r„/2)  cos(k  •  rJ2) . 


(3.45a) 


Here  we  made  the  common  assumptions  that  0^  =  d.-  d  and  that  the  system  is  quasi  one-dimensional 
with  /r,  +  A:,  perpendicular  to  the  lattice  vector  [55,  56).  Equation  (3.45a)  is  known  as  the  coherent 
initial  condition,  which  contains  site  populations  (m  =  n)  as  well  as  intersite  coherences  (m  n).  It  is 
expected  to  be  most  relevant  in  the  case  of  resonant  pumping.  Wong  and  Kenkre  [93]  have  also 
introduced  the  diagonal  initial  condition 


(b:(0)5„(0))  =  |i?|^cosUtJ2)5„„, 


(3.45b) 


which  only  contains  populations.  This  condition  is  important  in  the  case  of  off-resonance  excitation 
followed  by  ultrafast  vibrational  relaxation  [4]. 

We  will  restrict  our  treatment  to  the  case  of  coherent  excitation.  Combining  eqs.  (3.43).  (3.44c). 
(3.24b).  and  using  appendix  B,  we  eventually  arrive  at  the  Laplace  transform  of  the  population 
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W{k,-k,.s)^2d,  m{s-  \{J{k,)-Jik^)]  +  r^y} 


( 1  -  i  S  Is  -  iAp  -  t,/.2)  +  Utp  +  4/2)  +  r  *  y| 


This  can  be  rewritten  as 


m4.  -  4„ .) = 2  mH- ^  r  -  ^  I ) 


(3.4ba) 


(3.46b) 


.y(  p)^J{p-  kJ2)  -J(p  +  y 2) ,  '  (3.46c) 

where  the  A:g-dependence  of  i(  p)  is  suppressed  for  compactness. 

One  interesting  exact  result  follows  directly  from  eq.  (3.46):  in  the  absence  of  e.xciton  dispersion 
[/(A)  =  constant],  or,  in  particular,  in  the  absence  of  interactions,  the  population  W{k^  -k.,t)  only 
decays  with  the  trivial  rate  y.  independently  of  the  scattering  strength  t.  The  TG  signal  thus  decays  with 
rate  2y  (eq.  3.41a).  The  physical  explanation  is  that  all  coherences  Q{k,  p.  t)  now  have  the  same 
frequency  (zero),  so  that  they  do  not  dephase  with  respect  to  each  other  and  keep  adding  up  coherently 
to  the  population  (eq.  3.24b)  at  all  times,  irrespective  of  the  scattering  rates  between  them.  This  result 
is,  of  course,  not  restricted  to  the  Haken-Strobl  model. 

In  the  case  of  general  interactions,  it  is  hard  to  transform  eq.  (3.46b)  back  to  the  time  domain  and 
we  will  therefore  restrict  ourselves  from  now  on  to  a  study  of  the  long-time  behavior  of  the  grating 
signal.  Analytical  results  valid  for  all  times,  have  been  derived  by  Garrity  and  Skinner  for  a 
one-dimensional  system  with  nearest-neighbor  interactions  [55].  Their  important  conclusions  concerning 
the  limits  of  coherent  and  incoherent  motion  on  the  experimental  length  scale,  however,  can  also  be 
reached  from  our  expression.  We  first  note  that  the  long-time  limit.  \r+  y-  \J{  p)\t  >  1,  is  governed 
by  the  small  5  region  (|5|  1/’+  y  -  \J’{p)\)  in  the  Laplace  domain.  Equation  (3.46b)  then  yields 

t  ■<r«  >['^(  P)  ~  Pp)l 

W(k,-k..i)  =  2d,d*e\p{at).  a  =  i^(pJ-y  +  ^  S - —■  (3.47a.  b) 

iy  p  /  +  y  -  i.y(  p) 

In  the  remainder  of  this  section,  we  will  assume  that  -^{pj  =  J{k^)  -  /(it,)  =0.  which  holds  under 
very  general  symmetry  conditions  with  respect  to  the  experimental  set-up.  We  further  note  that 
■H  p)~  -k^  -  v^(  p).  with  v^ip)  the  exciton  group  velocity  defined  under  eq.  (3.38a).  Let  us  first 
concentrate  on  the  limit  of  weak  scattering:  J'{p)>  t  +  y  for  most  p  in  the  first  Brillouin  zone.  This  is 
the  coherent  limit,  as  the  exciton  scattering  length  v^{p)/r  h  now  much  larger  than  the  experimental 
length  scale  \kj[~\  We  then  obtain 

a  =  ~y  -  r .  (3.48a) 

so  that,  using  eq.  (3.41a),  we  find  for  the  intensity  of  the  TG  signal  normalized  to  its  initial  value 
5(t)  =  exp[-2(y  +  Or] .  (3.48b) 

This  result  has  also  been  found  by  Garrity  and  Skinner  [55].  In  the  extreme  case  of  /'  =  0.  no  grating 
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decay  is  observed,  apart  from  the  trivial  populafion  decay  exp(-2'yf).  This  result  is  rigorous  for  all 
times,  as  is  directly  seen  from  eq.  (3.46b),  and  is  easily  understood:  for  r  =  ()  the  coherent  initial 
condition  is  an  exact  eigenstate  of  the  system.  By  contrast,  the  diagonal  initial  condition  is  a 
superposition  of  eigenstates  and  leads,  even  in  the  case  T  =  0,  to  a  nontrivial  decay  of  the  grating  on  a 
time  scale  {k^v^y^  with  a  typical  exciton  group  velocity  (ballistic  motion)  [55,  93]. 

We  now  turn  to  the  opposite  limit  of  strong  scattering,  r>  y  and  r>  J'ip)  (all  p),  where  the  exciton 
scattering  length  is  much  smaller  than  the  grating  length  scale  {diffusive  or  incoherent  motion). 
Equation  (3.47b)  then  reduces  to 

a  = -y- DJ*,  -  ,  (3.49a) 

with  the  exciton  diffusion  constant  as  defined  in  eq.  (3.36).  In  the  derivation  of  eq.  (3.49a),  we  used 
that  Sp,^(p)  =  0,  as  is  explained  in  appendix  C.  We  now  obtain  for  the  signal  intensity  [55] 

5(t)  =  exp[-2(y  +  DjA,  -  k.\-)T] .  (3.49b) 

We  can  get  an  idea  how  well  this  long-time  expression  describes  the  actual  decay  by  using  in  eq.  (3.46b) 

y  _ 1  y  f 

P  s-i^{p)  + r+y~'  (5+r+y)' 7^  ■ 

(Equation  (3.49)  agrees  with  taking  5  =  0  in  the  r.h.s.  of  this  expression.)  The  thus  obtained 
approximation  to  W{k^  ~  5)  can  be  Laplace-inverted  analytically  and  leads  to  an  additive  correction 

in  eq.  (3.49b).  This  correction  term  has  a  relative  magnitude  which  is  much  smaller  than 

unity,  and  decays  on  a  time  scale  t  \  which  is  very  fast  compared  to  the  decay  in  eq.  (3.49b).  This 
suggests  that  eq.  (3.49b)  is  a  good  approximation  to  the  actual  signal  over  the  entire  observable  decay. 
Furthermore,  we  note  that  the  present  result  is  not  affected  by  the  exact  initial  condition  and  is  also 
found  for  diagonal  excitation  of  the  system  [55,  93].  The  reason  is  that  in  the  incoherent  limit,  intersite 
coherences  anyhow  relax  very  fast  compared  to  the  grating  decay  [cf.  eq.  (3.28b)]. 

We  will  now  discuss  the  result  eq.  (3.49)  in  more  detaiF.  The  signal  decay  rate  consists  of  a  trivial 
contribution  due  to  population  relaxation  (2y)  and  a  contribution  from  the  exciton  motion,  which  is 
proportional  to  |/c,  -  k,\'.  This  is  characteristic  for  diffusive  motion  [1.  2a,  55.  93]  and  leads,  for  small 
cross  angles  0  between  the  two  pump  pulses  (fig.  1).  to  a  linear  relation  between  the  observed  decay 
rate  and  0".  This  relation  allows  us  to  distinguish  experimentally  between  diffusive  and  coherent 
exciton  motion  on  the  scale  |fc,  -  ‘.  In  fig.  4  we  display  the  TG  signal  showing  incoherent  (diffusive) 

motion  and  its  dependence  on  the  fringe  spacing  through  the  variation  of  0.  In  fig.  4a,  the  TG  probes 
exciton  motion  in  disodium  fluorescein  in  ethanol  [8],  whereas  in  fig.  4b  it  probes  carrier  motion  in  a 
semiconductor  (GaAs).  In  both  cases  the  TG  decay  becomes  slower  as  0  is  reduced.  In  fig.  5  we  show 
the  0"  variation  of  the  TG  decay  rate  in  anthracene  at  two  temperatures  ( 10  K  and  20  K)  [2b].  Note  the 
diffusive  character,  even  at  these  low  temperatures.  In  spite  of  an  active  search,  coherent  motion  has 
never  been  directly  observed  in  transient  grating  experiments. 

In  the  incoherent  limit,  an  interesting  relation  exists  between  the  D4WM  and  TG  signals.  Namely, 
the  amplitude  of  the  D4WM  signal  (eq.  C.4)  can  be  obtained  by  evaluating  the  Fourier  transform  of  the 
amplitude  of  the  TG  signal  [exp(aT)]  at  the  frequency  w,,.  This  single-Fourier-transform  relation 


0.0  0.5  1.0  1.5  2.0  2.5  3.0 


Time  (ns) 

Fig.  4.  (a)  Fringe  spacing  dependence  of  the  TG  signal  probing  excited-state  energy  transport  in  disodium  fluorescein  in  ethanol,  for  (4/2  =  8°  and 
22”  (excitation  wavelength  266  nm)  (8|.  For  the  smaller  angle,  the  signal  is  dominated  by  acoustic  effects  (note:  time  scale  is  different  for  (4/2  =  8”). 
At  (4/2  =  22°.  the  acoustic-signal  contribution  damps  out  due  to  acoustic  attenuation,  (b)  Fringe  spacing  dependence  of  the  TG  signal  probing 
carrier  density  diffusion  in  GaAs.  =  2  x  10''*  photocarriers /cm’.  The  measured  ambipolar  diffusion  constant  is  in  good  agreement  with  the 
known  value  for  high  purity  GaAs  (0„  =  Scm’/sec).  Curve  1:  ©  =  4.9°,  curve  2;  ©  =  3°.  curve  3:  ©=1.1°  (see  ref.  1117|). 
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Fig.  5.  (a)  The  decay  rate  of  the  transient  grating  signal  versus  O'  for  two  anthracene  crystals  at  10  K.  along  the  a  axis  I’b).  The  magnitude  of  the 
slope  equals  Sir'DJA'.  with  A  the  wavelength  of  the  pump  pulses,  and  thus  yields  directly  the  exciton  diffusion  constant  D^.  The  value  of  the 
intercept  equals  2y,  where  y  is  the  population  relaxation  rate,  (b)  The  decay  rate  of  the  transient  grating  signal  for  the  same  two  crystals  at  20  K. 
The  difference  in  slopes  is  due  to  the  temperature  dependence  of  the  exciton  diffusion  constant.  As  the  temperature  increases,  a  decrease  in  the 
diffusion  constant  is  observed.  The  average  diffusion  constant  obtained  from  these  measurements  is  roughly  10  times  larger  than  the  expected  value 
for  incoherent  exciton  diffusion  (51|. 

between  nonlinear  frequency-domain  and  time-resolved  techniques  is  not  trivial  and  was  first  estab¬ 
lished  by  Loring  and  Mukamel  [36].  In  concluding  this  section,  we  mention  that  the  TG  signal  for  cases 
intermediate  between  the  coherent  and  incoherent  limits  has  been  studied  by  Garrity  and  Skinner 
(using  the  Haken-Strobl  model)  [55]  and  by  Wong  and  Kenkre  (using  a  generalized  master  equation 
approach)  [93].  Finally,  we  remark  that  we  only  studied  the  TG  within  the  two-particle  description.  The 
single-particle  (local-field)  description  would  lead  to  the  exponential  decay  eq.  (3.48b),  independently 
of  the  magnitude  of  t.  Diffusive  exciton  transfer  cannot  be  incorporated  at  this  level  of  description, 
which  is  intimately  related  to  the  fact  that  in  the  single-particle  factorization  no  resonance  with  respect 
to  is  found  in  the  D4WM  signal.  This  illustrates  the  limitations  of  the  local-field  approach. 


3.5.  The  Boltzmann  and  the  diffusion  equations  for  exciton  transport 

Further  physical  insight  into  the  equation  of  motion  (3.25)  for  the  two-particle  coherences  Q{k,  p.  t) 
can  be  obtained  by  making  the  connection  to  macroscopic  transport  equations,  such  as  the  Boltzmann 
equation  and  thr  diffusion  equation.  To  this  end,  we  define  the  Wigner  phase-space  distribution 
function  d>(r.  p,t)  by  [94.  95] 

S  Q{k,  p,  t)  exp(i*  •  r) . 


(3.50) 
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Here,  r  and  p  play  the  role  of  the  classical  position  and  momentum,  respectively,  of  the  particle  that  is 
described  by  the  Wigner  distribution.  Comparing  with  eq.  (3.24a).  it  is  clear  that  r  is  the  position  of  the 
exciton  ■  center-of-mass".  That  p  indeed  corresponds  to  the  exciton  momentum,  becomes  evident  if  we 
consider  the  pure  state  for  a  completely  delocalized  exciton  with  momentum  p,,;  then,  namely. 
Q{k,  p.i)~  The  Wigner  distribution  function  carries  the  full  quantum  behavior,  yet  it  has  a  well 

dehned  classical  analogue:  as  ^-^0  it  reduces  to  the  classical  density  matrix  of  the  particle  (exciton  in 
our  case),  and  it  is  therefore  very  useful  in  providing  a  semiclassical  insight  in  quantum  dynamics  [94]. 

In  this  section  we  will  only  be  interested  in  the  homogeneous  equation  of  motion  for  the  two-particle 
variables,  i.e..  we  ignore  the  last  term  in  eq.  (3.25).  which  couples  to  the  electric  held.  .Also,  we  will 
only  consider  a  purely  imaginary  self-energy  matrix  l{k:  p.  p  )  =  -iTfit;  p.  p  ).  The  Wigner  distribu¬ 
tion  then  obeys  the  equation 


^  (^(r.  p.  /  )  =  2  2!  J{a)  sin(  p-a)6{r^  ail.  p.  i)  - I'ia-.  p.  p')6(r  ^  a  1.  p’.t) . 


V  V 


f(a;  p,  p')=]j^  rik:  p.  p')  exp(-i/c  •  a.'2) 

k 


(3.51a) 

(3.51b) 


The  first  term  in  eq.  (3.5 la)  describes  coherent  exciton  motion  on  the  lattice  [we  used  J{r)  =  /(-r)];  the 
second  term  is  due  to  (phonon)  scattering.  For  scattering  kernels  that  do  not  depend  on  k, 
r{k:  p.  p')  =  g(p.  p').  we  have 


I'ia:  p.  p')  =  S^  „g(p.  p  ) .  (3.52) 

This  implies  that  the  position  is  not  affected,  whereas  the  momentum  is  scattered  as  in  a  discretized 
Boltzmann  equation  with  collision  kernel  g(  p.  p  ).  An  interesting  class  of  collision  kernels  is  found  by 
further  restriction  to 


g(p.  p')  =  (t  +  -  tg(p) .  (3.53) 

with  g(p)  =  1  (sum  over  the  first  Brillouin  zone).  The  Haken-Strobl  model,  cq.  (3.28).  is  contained 
within  this  class,  with  g(p)=  1/:V.  For  a  general  g(p).  the  equation  of  motion  now  reads 

d  V 

-j-  (Mr.  p.  f )  =  2  2^  J(a)  sin(  p-  a)(h(r  ^  a  1.  p.  t) 

-  I  'L  [k’(  P  )Mr.  p.n  -  ,e(  p)(b(r.  p  .  Mj  -  yMr.  p.i)  ■  (3.54) 

P 

The  second  term  in  this  equation  has  the  form  of  the  BGK  strong  collision  operator  in  the  Boltzmann 
equation  [96.  97|.  in  which  collisions  occur  with  rate  t  and  the  momentum  after  each  collision  is 
distributed  according  to  g(p).  In  this  model.  g(p)  is  the  equilibrium  momentum  distribution.  This 
implies  in  particular  that  the  Haken-Strobl  model  is  a  high-temperature  model,  as  the  equilibrium 
distribution  is  then  uniform  over  all  momenta.  The  strong  collision  operator  conserves  the  number  of 
particles  (population);  a  population  loss  with  rate  y  is  described  by  the  last  term  in  eq.  (3.54).  This  is 
also  nicely  illustrated  in  yet  another  representation  for  the  two-particle  variables,  which  is  found  by 
Fourier  transforming  d){r.  p.  t)  with  respect  to  p  [95] . 
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6(r.  s.  t)  = 


1 

^  ZMr.  p.  i)exp{-ip‘s) . 


(3.55) 


Comparison  with  eq.  (3.24a)  shows  that  s  plays  the  role  of  a  relative  coordinate  and  <6(r,  s  =  0.  t) 
represents  the  exciton  population  at  site  r. 

Equation  (3.54)  now  translates  into  |95], 

^  i{r.  s,  /)  =  -i^  y(fl)[<^(r+  all,  s  -  a,  t)  -  0(r  +  ail,  s  +  a,  r)] 

-  {r+y)Mr,  s,  t)  +  fg(s)4(r,0.  r) .  (3.56a) 

^(5)  =  Z^(p)exp(-ip-s).  ■  (3.56b) 

p 

As  a  consequence  of  the  normalization  of  g{p),  we  have  g(s  =  0)  =  1.  so  that  the  last  two  terms  of  eq. 
(3.56a)  indeed  describe  a  net  loss  of  population  with  rate  y. 

We  return  to  the  equation  of  motion  (3.54)  for  the  Wigner  distribution  and  direct  our  attention  to 
the  coherent  term.  We  will  assume  that  the  interaction  is  symmetric  {J{r)  =  /(-r)]  and  has  a  short  range 
R.  Furthermore,  we  assume  that  g{p)  is  centered  at  optical  momenta  p,  for  which  1.  In  the 

continuum  approximation,  eq.  (3.54)  then  becomes  [95] 

-  d{r.  p,  t)=  p,t)-  t  j  dp' [g{p')<b(r,  p.  t)  ~  g(p)d)(r.  p’.  /)]-  y<A(r.  p,  t) . 

(3.57a) 

where,  for  an  isotropic  d-dimensional  s,ystem, 

(m*)-'  (3.57b) 


Here,  n  runs  over  the  original  lattice  and  =  lr„|.  Furthermore,  gip)  is  defined  in  analogy  with  g{p), 
but  has  continuum  normalization  J  g(p)  dp  =  1.  Equation  (3.57a)  clearly  has  the  form  of  the  Boltzmann 
equation  for  a  classical  particle  of  mass  m*.  It  can  now  be  shown  quite  generally  that  in  the  limit  of  high 
friction,  t.  eq.  (3.57a)  reduces  to  a  diffusion  equation  for  the  particle  position  distribution  function 
defined  bv 


6(r.t)  =  ^  dpMr.  p.t) .  (3.58) 

Explicitly,  one  finds  (see  ref.  [97],  chapter  10  and  appendix  2): 

(dldt)d{r.  t)  =  DV'd)(r,  t)  -  y<b{r,  t) ,  (3.59a) 

with  the  diffusion  constant  defined  by 

D  -  i  (( ,,/m*  =  i  (  i  2  .  (3.59b) 

where  (•  •  •)^-  takes  the  average  over  the  equilibrium  distribution  g{p). 
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The  direct  solution  of  the  equation  of  motion  for  Q{k.  p,  t)  within  the  Haken-Strobl  model  also 
leads  to  a  diffusion  constant  (eq.  3.37), 

D,=iVdn'2Ar.)rl.  (3.60) 

n 

It  is  clear  that  in  general  D  does  not  coincide  with  even  if  we  take  g{p)  -  {allirf,  which  is  the 
equilibrium  distribution  in  the  Haken-Strobl  model  {a  is  the  lattice  spacing).  The  quantitative 
difference  between  the  two  above  diffusion  constants  can  be  illustrated  by  giving  their  values  in  the  case 
of  nearest-neighbor  interactions  (7)  on  a  lattice  with  spacing  a. 

D^  =  {2/r)J-a- ,  D  =  {47r-d/3r)J'a- .  (3.61a,  b) 

The  solution  to  this  apparent  discrepancy  is 'that  for  the  Haken-Strobl  model,  strictly  speaking,  the 
Boltzmann  equation  cannot  be  derived  through  the  above  method,  because  the  equilibrium  distribution 
g{p)  is  uniform  instead  of  centered  at  optical  wave  vectors.  Consequently,  eq.  (3.59b)  cannot  be 
derived  within  the  Haken-Strobl  model.  Alternatively,  for  a  very  broad  equilibrium  distribution  g(p), 
the  classical  picture  of  a  momentum-independent  effective  mass  m*  is  an  oversimplification. 

That  it  is  possible,  however,  to  derive  in  a  rigorous  way  a  diffusion  equation  for  the  exciton  motion 
within  the  incoherent  limit  of  the  Haken-Strobl  model  is  clear  from  the  discussion  following  eq.  (C.5). 
A  very  elegant  derivation  has  been  put  forward  by  Reineker  and  Kiihne  [54,  98.  99).  These  authors  first 
derive  the  Pauli  master  equation  for  the  site  exciton  populations  in  the  strong  scattering  limit.  In  our 
notation,  this  equation  reads  [54.  98[ 

dW 

H..,  -2|y(r,„)|’/r.  (3.62a,b) 

Ui 

Equation  (3.62a)  has  been  extensively  studied  within  the  field  of  incoherent  energy  transfer  [100]; 
is  the  well-known  Forster  rate  of  energy  transfer  between  the  molecules  n  and  n'  [41.  101).  For  dipolar 
interactions,  has  the  characteristic  dependence.  In  the  continuum  approximation,  the 

diffusion  equation  for  the  exciton  population  density  is  now  easily  derived  from  eq.  (3.62a).  by 
expanding  around  the  lattice  point  r„  [54,  99).  For  an  isotropic  system  the  diffusion  constant 
obtained  in  this  way  indeed  coincides  with  eq.  (3.60). 


4.  Nonlinear  optical  response  of  polaritons 

4.1.  Canonical  transformation  and  operator  equations 

In  this  section  we  develop  a  polariton  theory  of  nonlinear  optical  response  for  the  case  that  the 
radiation  and  the  matter  degrees  of  freedom  are  strongly  coupled  and  the  TG  and  D4WM  experiments 
are  dominated  by  polariton  dynamics.  Polaritons  are  the  combined  radiation-electronic  eigenmodes  of 
the  crystal,  which  have  a  dispersion  diagram  that  can  differ  profoundly  from  that  for  excitons  (see  fig.  6) 
[57,  58].  Of  course,  the  optical  response  of  a  crystal  is  in  general  determined  by  its  proper  elementary 
excitations,  so  that  optical  signals  must  exhibit  resonances  and  broadenings  determined  by  the 
energetics  and  dynamics  of  polaritons  instead  of  excitons.  The  theory  of  section  3  cannot  account  for 
those  effects:  the  calculation  of  response  functions  and  susceptibilities  was  based  on  a  factorization  of 
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the  electric  field,  thus  treating  it  effectively  as'an  external  c-number  and  neglecting  part  of  the 
correlations  between  the  material  (excitons)  and  the  radiation  fields.  Such  a  scheme  leads  to  signals 
characterized  by  the  exciton  dispersion  and  dynamics.  Although  it  is  in  principle  possible  to  define  and 
calculate  susceptibilities  without  factorization,  this  does  not  lead  to  a  very  practical  scheme  [66.  67], 
Here,  we  rather  follow  an  approach  in  which  a  hierarchy  of  nonlinear  equations  of  motion  for  polariton 
variables  inside  the  crystal  is  built,  which  directly  leads  to  the  optical  signal  (polariton  hierarchy)  [76]. 
This  method  is  described  in  more  detail  in  the  next  section. 

As  in  the  exciton  theory,  our  approach  starts  by  deriving  a  basic  set  of  operator  equations.  Now, 
however,  we  also  consider  equations  for  the  photon  operators,  as  the  radiation  field  is  treated  as  an 
explicit  degree  of  freedom.  We  use  eq.  (2.9)  as  starting  point  and  neglect  Umklapp  contributions  to  £*. 
Furthermore,  we  note  that,  at  every  wave  vector  k.  we  only  need  to  consider  photons  with  polarization 
in  the  plane  spanned  by  k  and  the  molecular  transition  dipole  ft  :  photons  polarized  perpendicular  to 
this  plane  do  not  couple  to  the  crystal  and  form  no  new  eigenmodes  with  the  excitons.  We.  thus,  drop 
the  polarization  labels  A  on  the  photon  creation  and  annihilation  operators  and  arrive  at  (all  operators 
taken  at  time  t) 


-iQ 

(I 

-iQ 

1  ‘h  ] 

1  d 

1  iC, 

-n-D^  +  i/U) 

-iQ 

Bk 

i  d/ 

0 

1 

-iQ 

kc 

-iQ. 

i  1 

“  ki 

0, 

iQ 

n  -  D,  -  i 

mj 

Mk) 


/M 

1 

0 

-1/ 


(4.1a) 

Inp-- 

A(^)=  -c,:^).W(k-k')].  .  (4.1b) 

-  * 

C\  =  (Ivkcph  'p.'  )'  ■  .  (4.2a) 


Dj  =  J(k)  +  'inpfi  'p'  . 


p  =  p  - ip • kf  k  . 


(4.2b.  c) 


In  the  derivation  of  eq.  (4.1)  we  used  eqs.  (2.6)-(2.8)  to  express  E  in  terms  of  photon  and  exciton 
creation  and  annihilation  operators.  Of  course,  the  equations  for  and  B  ^  arc  equivalent  to  eq. 
f.T  1 1.  The  inclusion  of  l'{k).  the  phonon-induced  exciton  damping,  is  not  fully  consistent  here,  as  it  was 
calculated  without  accounting  for  polaritons.  We  will  come  back  to  this  further  on  in  this  section.  From 
eq.  (4.1)  one  easily  derives  the  Maxwell  equations  in  the  electric  dipole  approximation  [71.  76] 


E  {k.i)=  -  -  -  A  {k,t). 

C  ()l 


k-  ^  —-]E-(k.()=^  -  ~  ^  P  ik.i) . 

(■'  ctr  c 


(4.3a.  b) 


These  equations  are  not  affected  by  the  nonlinear  term  -S  {k)  in  eq.  (4.1). 

We  now  consider  eq.  (4.1)  in  the  Bose  approximation  X{k)  =  i).  It  is  then  a  closed  4x4  matrix 
equation  that  completely  describes  the  linear  propagation  of  electromagnetic  waves.  The  eigenmodes  of 
this  proolem  are  the  polaritons.  which  are  related  to  the  photons  and  excitons  by  a  canonical 
transformation  [57.  58.  61 1.  We  follow  Hopfield's  notation  and  search  for  a  transformation  [57] 
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such  that 


(4.4) 


(d/d/)^4,  = -ioj^  ( Bose  approximation! .  (4.5) 

Here  denotes  the  annihilation  operator  tor  a  polariton  with  wave  vector  k  in  branch  /■.  and  is  the 
frequency  ot  this  polariton.  In  combination  with  eq.  (4.1)  this  defines  an  eigenvalue  problem  whose 
secular  equation  gives  the  polariton  dispension  relation. 


iUr 


to 


ki 


=  1  + 


JJ _ 

i/U)l'  ^  /i;  ■ 


/  =  (iSTTp/i"  ii.  ft)'  ■  . 


(4.6a.  b) 


where  is  the  exciton  dispersion  eq.  (.v6).  .As  the  polariton  dispersion  is  by  definition  equivalent  to 
the  one  for  electromagnetic  waves  in  the  crystal,  (kc  must  also  define  the  transverse  dielectric 
‘unction  t  (A:.(a*^,).  Comparison  with  cq.  (3.12b)  shows  that  the  thus  obtained  transverse  dielectric 
function  is  identical  to  the  one  found  through  exciton  response.  The  linear  response  ot  the  crystal,  apart 
from  a  possible  chani’e  of  the  damping  constant  /  (A)  |sw  helon  ],  is  not  atfccted  hv  the  introduction  of 
the  polariton  koncept  [57.  58.  SOj. 

We  now  first  analyze  the  polariton  dispersion  and  transformation  in  more  detail  for  the  important 
case  /’(A)  =  ().  i.e..  totally  neglecting  the  phonon  bath.  This  is  the  standard  case  considered  in  the 
literature.  Equation  (4.6)  then  gives  the  usual  dispersion  diagram  (cf.  fig.  h)  with  two  branches 
separated  by  the  polariton  stopgap:  the  frequency  interval  where  no  real  wave-vector  solutions  to  the 
dispersion  relation  exist,  so  that  no  waves  with  those  irequencies  can  propagate  m  the  crystal.  In  an 
atomic  crystal,  the  stopgap  ranges  from  the  transverse  («>  )  to  the  longitudinal  (w  )  exciton  frequency 
|43.  .■'7.  58|.  For  our  crystal  of  two-level  molecules,  these  boundaries  depend  on  the  direction  of 


tig.  Upii-.il  P‘'ljrilon  dispcrvii'n  curves  in  the  optie.il  reeu'n  li>r  .in  .iti'mie  ervsuil  uhick  solid  curves).  I  he  di.ieon.ii  line  represents  the  pure 
photon  dispersion  curve  (w,  =  3c  I  The  shaded  region  between  the  transverse  \w  l  and  longitudinal  lui  I  crvsial  evciion  Irequencies  is  the  stopgap, 
where  no  polariton  modes  exist  ti  indicates  the  atomic  transition  trequencv.  For  a  crvsial  of  two-level  molecules,  the  stiipeap  position  depends  on 
the  direction  cif  propagation  (see  text) 
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propagation  and  are  replaced  by  and  ({};,  *  / "  respectively,  with  II,  =  linij^,,  and  /  as  defined 
in  eq.  (4.6).  (Exciton  dispersion  is  neglected  when  considering  the  stopgap.) 

The  strongest  mixture  between  excitons  and  photons  occurs  in  the  photon-  and  exciton-band  crossing 
region  kc  ^  =  Q.  In  that  region  the  frequency  difference  between  the  two  branches  reaches  its 

minimum  value,  which,  to  a  good  approximation,  is  given  by  /[60b].  From  eq.  (4.6).  we  observe  that  / 
is  a  measure  for  the  density  of  oscillator  strength  in  the  crystal  and  in  combination -with  eq.  (4.2a).  it  is 
also  seen  that  fl2  is  the  coupling  Q  between  the  exciton  and  the  radiation  field  in  the  band-crossing 
region.  For  molecular  crystals,  we  usually  have  f  <11.  For  example,  for  the  lowest  singlet  a-exciton  in 
naphthalene,  we  have/=45cm  '  and  Z},,  ^/2  =  31  500  cm  '’  [60b);  for  the  same  exciton  in  anthracene, 
we  find  / ~  1(K)0  cm  ‘  and  ==  25  0(X)  cm  ‘  [bOa).  ( Both  /  values  acco  mt  for  the  crystal's  background 
index  of  refraction.) 

The  polariton  transformation  coefficients,  which  are  determined  by  the  eigenvectors  of  the  diagonali- 
zation  procedure,  also  nicely  demonstrate  the  mixed  character  of  the  polaritons.  In  the  case  Fik)  =  0. 
the  coefficients  are  found  to  be 
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(a»  -f-  kc){oj'  -  Ql  ~  f') 
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(4.7a) 


(4.7b) 


(4.7c.  d) 


where  ip  is  an  arbitrary  phase  and  oi  stands  for  oij, .  These  coefficients  have  been  normalized  according 
to  [57.  58.  61] 


In  spite  of  the  complexity  of  eqs.  (4.7),  several  general  and  instructive  conclusions  can  be  drawn. 
First,  the  coefficients  and  usually  have  absolute  values  small  compared  to  unity  (see.  e.g.,  ref. 
[68]).  Second,  the  limits  for  small  and  large  wave  vectors  are  easily  studied.  In  the  limit  iA:|— »0.  we  find 
for  the  upper  branch  |.r|  =  1,  |  vv|  =  \y\  =  |r|  =  0  (up  to  order  fill),  so  that  the  polariton  is  a  pure  exciton 
there.  For  |itl— the  same  branch  represents  a  pure  photon  (|wi  =  1,  |.v|  =  l.vi  =  i:|  =6).  This  could, 
of  course,  have  been  already  guessed  from  the  dispersion  diagram  (fig.  6).  Equally  natural  results  are 
found  for  the  lower  branch.  For  intermediate  wave  vectors  the  polaritons  smoothly  change  character 
from  exciton-like  to  photon-like  (or  vice  versa).  Third,  at  the  band  crossing  (kc  —  II).  we  find 
Kl"  “  kl’  -  ;  for  both  branches,  confirming  that  the  strongest  photon-exciton  mixture  occurs  there. 

We  note  that  the  coefficients  (4.7)  do  not  e.xactly  coincide  with  the  ones  obtained  within  the  minimal 
coupling  Hamiltonian  [57.  58).  but  their  basic  features  are  not  affected  by  this.  The  dispersion  relation 
(4.6).  which  is  an  ‘observable",  is  of  course  independent  of  the  choice  of  the  Hamiltonian.  The  present 
polariton  transformation  does  coincide  with  the  one  derived  in  ref.  [102].  if  in  the  latter  we  ignore  the 
coefficients  coupling  to  the  higher  Brillouin  zones  (Umklapp  contributions)  and  neglect  the  retardation 
in  the  effective  interatomic  interaction.  Finally,  we  will  need  the  inverse  polariton  transformation, 
which  can.  quite  generally,  be  shown  to  read  [61] 


We  now  return  to  the  general  polariton  transformation  with  r(k)  ¥=  0.  The  two  branches  found  above 
are  then  no  longer  the  correct  eigenmodes.  For  //T(A:)<g  1.  the  maximum  amount  of  mixing  between 
photons  and  excitons  is  of  the  order /7r(A:),  so  that  the  polaritons  are  to  a  good  approximation  just  the 
photons  and  (damped) -excitons  [60.  68].  In  other  words,  if  the  exciton-radiation  coupling,  /.  is  small 
compared  to  the  coupling  of  excitons  to  phonons  (or  any  other  perturbing  degree  of  freedom  [103]). 
measured  by  r{k).  no  strongly  mixed  radiation-matter  eigenmodes  exist.  We  then  expect  the  exciton 
description  of  optical  response  to  be  accurate.  We  note  that  in  this  case  of  strong  exciton  damping,  the 
elimination  of  the  phonon  bath  [leading  to  Fik)]  before  considering  the  coupled  exciton-photon  svstem 
is  fully  justified. 

Conversely,  the  condition  f/r{k)>  1.  which  is  typically  met  in  low-temperature  pure  crystals,  marks 
the  region  of  strong  polariton  effects  in  optical  response.  It  is  then  inconsistent  to  account  for  phonons 
(or  other  bath  variables)  by  including  the  exciton  damping  r{k)  in  eq.  (4.1).  Instead,  we  should  first 
diagonalize  the  bare  exciton-photon  system  [T(^)  =  0]  and  then  account  for  other  degrees  of  freedom 
by  considering  their  perturbation  to  the  thus  found  polariton  eigenmodes  (polariton  self-energies)  [60b. 
68].  In  the  remainder  of  this  review  we  will  assume  that  the  limit  of  strong  polariton  effects  applies,  so 
that  from  now  on  the  polariton  dispersion  and  transformation  coefficients  refer  to  the  case  r{k)  =  0. 
The  polariton  self-energies  are  addressed  in  appendix  A  and  incorporated  in  the  equations  of  motion  in 
section  4.2.  For  completeness  we  note  that  the  theoretical  treatment  of  the  intermediate  region  where 
the  exciton-photon  coupling  roughly  equals  the  exciton-phonon  coupling,  is  very  complicated,  because 
it  is  then  impossible  to  indicate  a  simple  set  of  eigenmodes  that  are  weakly  perturbed  by  the  phonon 
bath. 

So  far  we  only  considered  polaritons  in  the  Bose  approximation,  which  results  in  a  purely  linear 
optical  response.  However,  using  the  above-obtained  polariton  transformation,  the  full  equation  of 
motion  (4.1)  yields 

d^*./dt=  +  i(x*,  -  z,jMk) .  (4.10) 

and  its  hermitian  conjugate  for  the  polariton  creation  operator  .  Equation  (4.10)  is.  apart  from  the 
omission  of  the  phonon-induced  self-energy,  exact  and  fully  accounts  for  nonlinear  behavior.  The  linear 
part  of  the  equation  is.  of  course,  equivalent  to  eq.  (4.5);  the  nonlinear  part  is  determined  by  eq.  (4.  lb) 
and  may  be  translated  into  polariton  creation  and  annihilation  operators  by  using  the  inverse  polariton 
transformation  (4.9).  A  typical  term  in  S'{k)  reads 

S  S  <.(«'*■. „  ,4J|| 

k  k  i'  »■  I-  ' 


'  ith  a(kk'k".  a  complicated  function  of  wave  vectors  and  branches.  Finally,  as  the  exciton 

population  operator  plays  a  crucial  role  in  the  nonlinearity,  it  is  useful  to  give  its  full  expression  in  terms 
of  polariton  operators.  Using  eqs.  (3.5)  and  (4.9)  we  arrive  at 
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4.2.  The  polariton  hierarchy  in  the  two-particle  description 

Before  turning  to  the  treatment  of  the  TG  and  the  D4WM  experiments  in  the  limit  of  strong 
polariton  effects,  we  first  discuss  the  general  approach,  the  polariton  hierarchy  [76].  Consider  a  finite 
crystal  that  is  still  large  enough,  however,  so  that  its  eigenmodes  are  to  a  good  approximation  the 
polaritons  discussed  in  section  4.1.  The  nonlinear  evolution  of  the  isolated  polariton  system  is  in 
principle  fully  described  by  eq.  (4.10),  but  thfe  equation  does  not  yet  describe  the  generation  of 
polaritons  by  the  external  laser  fields.  It  is  important  that  this  generation  is  not  caused  by  a  direct 
interaction  term  between  electromagnetic  fields  and  polaritons  in  the  Hamiltonian  [69],  because  the 
electromagnetic  field  itself  is,  inside  the  crystal,  fully  contained  within  the  polaritons.  Instead,  the 
coupling  occurs  through  the  boundary  conditions  at  the  surface  of  the  crystal.  If  an  external  electric 
field  with  wave  vector  /c'  and  frequency  {cj-  =  jifc'Jc)  is  incident  on  the  crystal,  it  launches  a  polariton 
which  is  specified  by  matching  the  boundary  conditions  for  the  expectation  values  of  its  electromagnetic 
field  components  to  the  external  field  (accounting  for  reflection).  Thus,  the  polariton  will  also  have 
frequency  its  wave  vector  kj  is  determined  by  the  component  parallel  to  the  crystal’s  surface 
(*y.il  =  4:'  ||)  and  by  its  magnitude,  [4:^1,  which  is  fixed  through  and  the  polariton  dispersion  relation 
(i.e.  the  crystal’s  dielectric  function)  [69].  This  agrees  with  Snell’s  law.  Finally,  the  amplitude  of  the 
polariton,  (ltj,(/)),  is  (in  the  frame  rotating  with  frequency  ta^)  proportional  to  the  external  field 
amplitude  £“. 

The  complete  solution  to  the  boundary  value  problem  is.  in  particular  if  spatial  dispersion  is  taken 
into  account,  a  complicated  problem  [43.  104-107].  Instead  of  rigorously  solving  this  problem,  we  will 
incorporate  the  above  well  accepted  ideas  about  the  generation  of  polaritons  by  adding  simple  source 
terms  to  our  equations  of  motion  (see  below).  In  the  absence  of  the  nonlinear  term  in  eq.  (4.10),  the 
thus  generated  “first-order"  polaritons  constitute  the  only  response  of  the  crystal  to  exciting  fields. 
Clearly,  however,  the  nonlinear  term  causes  fusion  of  first-order  polaritons  to  higher-order  ones  at  new 
(not  fundamental)  wave  vectors  and  frequencies.  Such  a  higher-order  polariton  generates  an  observable 
nonlinear  signal  outside  the  crystal  which  is.  again,  found  by  matching  boundary  conditions.  In 
particular,  the  signal  amplitude  will  be  proportional  to  the  amplitude  (ft,.  (0)  of  the  higher-order 
polariton.  The  problem  of  calculating  the  optical  response  of  the  crystal  now  boils  down  to  calculating 
the  amplitude  (It  „ ) ,  which  can  be  done  using  a  truncated  hierarchy  of  equations  of  motion  for 
polariton  expectation  values  based  on  eq.  (4.10). 

To  apply  this  scheme  to  the  TG  and  D4WM  experiments,  we  use  the  truncation  that  also  proved 
useful  in  the  exciton  theory,  namely,  we  factorize  the  exciton  population  in  the  e.xpectation  value  of 
\'(4c)  (two-particle  description).  From  eq.  (4.10).  we  then  obtain  (t\„(f)  =  (lt,  (0)l 

d  V 

^  ^*.,(0  =  -  1  Z  (4:)^*,,(r) 

+  (a:*.,  -  2  [7(4:V)£,,.,  (f) -  y*(*V')£ !,,,„(r)]VF(4: -  k'.  t) 

k'  „■ 

+  S  1(4:1^.  it;  )£;(/)  exp(-ia;/) .  (4.13a) 

/ 

Here.  (4c)  is  the  phonon-induced  polariton  self-energy  given  in  eq.  (A.l4b)  and 
y{ku)  =  C,(wl  +  yl)  +  \D,{xl  -  z*J  . 


(4.13b) 
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Finally,  the  last  term  in  eq.  (4.13a)  describes  the  above  discussed  generation  of  polaritons.  where  £y(t) 
is  the  (slowly  varying)  amplitude  of  the  jth  external  electric  field.  The  proportionality  constant  A{ku,  k'j) 
is  nonzero  only  if  the  jth  external  field  mode  exactly  matches  the  polariton  ku:  its  magnitude  is  then  a 
complicated  function  of  geometry  and  frequency  which  we  do  not  calculate  here.  Of  course,  this 
“source  term”  only  makes  sense  for  pulses  which  are  long  compared  to  an  optical  cycle. 

As  next  step  in  the  hierarchy,  we  need  equations  for  the  two-particle  variables  contained  in 
{W{k  -  k\  /)).  In  section  3  we  considered  the  exciton  two-particle  variables.  Q(k.  p,  t)  [cf.  eqs.  (3.24)], 
but  this  choice  is  not  useful  here,  because  these  variables  couple  to  expectation  values  of  products  of 
exciton  and  radiation  operators,  which  we  cannot  factor  in  the  limit  of  strong  polariton  effects.  Instead, 
it  is  more  appropriate  and  more  direct  to  consider  polariton  two-particle  variables.  From  eq.  (4.12),  we 
see  that  four  types  of  such  variables  are  involved.  However,  as  in  general  |z*J  ^  1  and.  moreover,  for 
the  self-energy  models  to  which  we  will  confine  ourselves  later  on,  |2*^|  we  will  only  maintain 

the  first  r.h.s.  term  in  eq.  (4.12),  so  that 


Wik.  t)  -  (2//V)  I  S  (*.  P.  0 


(4.14) 


(4.15) 


The  next  step  in  the  hierarchy  now  results  in  the  equation  of  motion  for  the  variables  ^^y  {k,  p,  t)f 
which  reads 

(d/dr)5,„.(A.  p,  r)  =  -  Wp,*  :.,,)r,,„  (*.  p.  /)  -  i  S  2  -(^.  i'-  P-  p')H,,„  (*.  p'.  0 

p 

+  ^p-k:i.A0  2  T(p  +  kl2.  t/,  *')£"(/)  exp(-ia>^7) 


+  ^p-*  :..  (0  2  A*(  p  -  kl2.  p'.  k]  )£'■',  Ui)  exp(iw,0  . 


(4.16) 


The  first  term  in  this  equation  reflects  the  linear  coherent  motion  [first  r.h.s.  term  in  eq.  (4.10)].  The 
nonlinear  coherent  motion  is  neglected  on  this  level,  because  it  couples  to  variables  which  are  of  fourth 
and  higher  orders  in  the  external  field  amplitudes,  which  is  beyond  the  leading  order  (3)  for  four-wave 
mixing  processes.  The  second  term  in  eq.  4.16)  accounts  for  phonon  scattering.  The  calculation  of  the 
self-energy  kernel  l(k,v\p.p')  [p  =  (p.  p'.  p".  p"')\  is  addressed  in  appendix  .A  (eq.  A. 16),  and 
remarks  similar  to  those  following  eq.  (3.26)  for  the  exciton  self-energy  apply  here.  In  particular,  we 
see  again  that  the  “center-of-mass”  wave  vector  k  is  conserved  in  the  scattering  process.  Finally,  the  last 
two  contributions  in  eq.  (4.16)  account  for  creation  of  two-particle  coherences  from  existing  single¬ 
particle  amplitudes  and  new  polaritons  generated  bv  external  fields  [cf.  the  generation  term  in  eq. 
(4.13a)]. 


4.3.  Polariton  dynamics  probed  by  transient  grating  spectroscopy 

We  now  apply  the  two-particle  polariton  hierarchy  [eqs.  (4.13a)  and  (4.16)]  to  the  TG  experiment  in 
the  limit  ot  strong  polariton  effects.  Two  external  pump  fields  [(itj.  w, )  and  ik'..  oi, )]  create  polaritons, 
/c,  p^  and  k.p..  resulting  in  a  grating  characterized  by  the  wave  vector  it,  -  k^.  After  a  delay  time  t,  a 
probe  pulse  (*(,  (o^)  generates  probe  polaritons,  k^p^,  which  interact  nonlinearly  with  the  grating  and 
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give  rise  to  signal  polaritons  with  wave  vector  k,  ^k^  -  ifc,  +  kj.  The  amplitude  ^  (t)  determines  the 
amplitude  of  the  observable  field  at  k'^  =  k\-  k’.,  +  k'^.  From  eq.  (4.13)  we  observ'e’that  the  nonlinear 
source  term  for  the  signal  polaritons  to  third  order  in  the  external  field  amplitudes  is  given  by 


(4.17) 


where- superscripts' in  parentheses  indicate  the  order  in  powers  of  the  external  field  amplitudes.  We  will 
assume  an  off-resonance  (square)  probe  pulse  that  is  short  compared  to  the  typical  time  scale  for 
population  decay.  On  this  time  scale,  the  amplitude  ^*'|,j(t)  then  behaves  as  a  delta  function  centered 
around  r  =  t,  so  that  the  time-integrated  signal  intensity  (to  third  order  in  the  external  field  intensities) 
is  given  by 


S(t)  .  I  dr  I f i;ijr)p  « I  r)H . 


(4.18) 


Thus,  as  in  the  exciton  theory  of  section  3.4,  the  TG  signal  is  determined  by  the  exciton  population  at 
the  grating  wave  vector  k^  =  k^-  k^.  This  is  a  direct  consequence  of  the  fact  that  the  exciton  population 
operator  is  realty  the  only  nonlinearity  in  our  system  (cf.  section  3.1;  the  photons  are  perfect  bosons). 
The  difference  with  the  theory  of  section  3  is,  of  course,  that  now  the  evolution  of  the  exciton 
population  occurs  within  the  complete  polariton  system  [eqs.  (4.14)-(4.16)]  instead  of  the  isolated 
exciton  system.  We  thus  take  into  account  that,  even  though  there  are  no  external  fields  during  the 
pump-probe  delay,  the  (microscopic)  radiation  field  still  affects  the  material.  The  TG  signal  is  a  probe 
for  polariton  dynamics  during  the  delay  period.  A  second  difference  between  the  exciton  and  polariton 
approaches  is  the  initial  condition  immediately  following  the  pump  pulses.  In  section  3,  these  were 
exciton  amplitudes  and  coherences  (eqs.  3.44),  whereas  in  ‘.he  present  treatment  polariton  amplitudes 
and  coherences  are  appropriate.  Using  eq.  (4.13a)  for  square  excitation  pulses  with  amplitudes  and  a 
duration  that  is  short  enough  to  neglect  the  self-energy  (7;,l-2'„„  |  -^  1),  we  obtain 


^k^v^iO  =  Mkit'i,k’)tE)expi-iu)it),  y  =  l,2,  0<r<T, 


(4.19) 


This  leads  in  particular  to  the  initial  condition  for  the  delay  period 


exp(-ia;^Tj 


(4.20) 


More  directly  relevant  to  the  description  of  the  TG,  are  the  initial  conditions  for  the  two-particle 
variables  s[lXk^,p,t).  Using  eq.  (4.19)  in  eq.  (4.16)  and  neglecting  the  self-energy  contributions 
during  the  pulses,  we  obtain 


P,  Tj  =  A{k,  u,,k[)A*{k,u,,  k',)E]El*T;  exp(i(o>,  -  cu,)tJ  6^ 


(4.21) 


with  p^  =  {k^-\-  kj)!!.  Note  that  this  result  is  consistent  with  the  factorization  (0^*,.(0)^‘’ = 
(l*i-  (0)*'’(^*,.(0)*'’  during  the  pulses,  which  results  from  the  fact  that  phonon  scattering  can  be 
neglected  for  these  early  times.  An  analogous  factorization  has  been  used  in  section  3.4  to  determine 
the  initial  two-particle  exciton  variables. 
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In  summary,  the  TG  signal  can,  in  principle,  be  calculated  from  eqs.  (4.18)  and  (4.14)  after  solving 
the  coupled  equations  of  motion  (4.16)  for  the  two-particle  variables  5'I,^,'(A:^,,  p.  i)  (for  all  p  and  with 
k^  =  k^-  k^  fixed)  during  the  pump-probe  delay  period.  During  this  period,  the  last  two  terms  in  eq. 

(4.16)  are  absent  (no  external  fields)  and  the  initial  condition  for  this  period  is  given  by  eq.  (4.21).  In 
practice,  we  have  to  restrict  the  treatment  to  simple  models  for  the  self-energy  (or  scattering)  kernel  to 
find  analytic  expressions  for  the  signal.  Before  discussing  such  models,  we  give  a  general  result.  In  the 
absence  of  scattering  (low  temperature),  the  TG  signal  does  not  decay.  Namely,  the  solution  of  eq. 

(4.16)  for  k  =  k^  is  in  this  case  trivially  given  by  the  initial  condition  (4.21)  with  an  additional  phase 
factor  exp(i(w.  -  w,)t].  As  this  factor  gives  the  only  r  dependence  of  W(k^,  t).  the  signal  does  not 
depend  on  t  at  all.  It  should  be  noted  that  this  result  does  not  rely  on  the  approximation  eq.  (4.14).  but 
also  holds  if  the  complete  form  eq.  (4.12)  for  the  exciton  population  is  used.  The  physical  explanation 
is,  of  course,  that  the  initially  created  polariton  coherence  (eq.  4.21)  is  an  eigenstate  of  the  system  in 
the  absence  of  scattering. 

We  now  discuss  two  specific  models  for  the  phonon-induced  scattering  kernel  I(k.  v.  p.  p')  in  eq. 

(4.16) ,  for  which  analytic  results  can  be  obtained.  In  the  exciton  theory,  we  used  the  Haken-Strobl 

model.  As  discussed  in  section  3.5,  this  is  a  (high-temperature)  strong-collision  model  that  scatters  all 
excitons  into  each  other  with  equal  rate  t,  irrespective  of  their  energies.  Using  the  same  model  for 
polaritons  would  be  very  unrealistic.  First,  because  these  excitations  span  an  enormous  bandwidth. 
Second,  it  is  only  the  exciton  component  of  the  polariton  that  couples  directly  to  phonons,  so  that 
photon-like  polaritons  cannot  be  scattered  very  strongly.  A  simple  model  that  takes  these  considera¬ 
tions  into  account,  is  a  restricted  Haken-Strobl  model,  which  is  defined  as  follows.  All  polaritons  with 
exciton  components  |.r*J"  >  ;  scatter  into  each  other  with  equal  rate  whereas  all  other  (photon-like) 
polaritons  do  not  take  part  in  the  scattering  process  at  all.  This  model  boils  down  to  the  Haken-Strobl 
model  for  the  restricted  set  of  exciton-like  polaritons.  which  (section  4. 1 )  constitute  the  upper  polariton 
branch  for  \k\  <k,  and  the  lower  branch  for  |it|>A:,,  marks  the  photon-exciton  hand 

crossing).  We  thus  consider  only  one  effective  branch  (i.e..  only  one  polariton  per  wave  vector)  with 
width  /(fig.  7),  so  that  branch  labels  can  be  dropped  without  ambiguity. 


Fig.  7.  Illustration  ot  the  restricted  Haken-Strobl  scattering  model  for  polantons.  All  polaritons  within  the  shaded  region  ol  width  /  (eq.  4  o)  have 
exciton  components  greater  than  4  (section  4.1)  and  are  assumed  to  scatter  into  each  other  with  equal  rates.  Sl„  is  the  exciton  frequency  at  zero 
(optical)  wave  vector  and  k„  =  tijc  indicates  the  exciton- photon  (dashed  diagonal)  band  crossing.  The  polariton  stopgap  is  not  shown  and  exciton 
dispersion  has  been  neglected. 
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In  the  Laplace  domain,  eq.  (4.16)  for  k  =  now  takes  the  form 

P,  ln)l='^’(*,.  P.  s) 

+  (/l/AT)  2  (>'.  s)  +  ^  ( =  0) . 


with  the  initial  condition 


(4.22a) 


=>='(*,,  p.l  =  0)«E\El'S 


(4.22b) 


Here  we  assumed  that  initially  exciton-like  polaritons  are  excited.  In  eq.  (4.22a)  we  added  an  overall 
loss  of  polaritons  with  rate  y^.  The  equation  of  motion  is  identical  to  eq.  (3.43),  except  that  the  exciton 
dispersion  and  scattering  rates  are  replaced  by  those  for  the  polaritons.  The  solution  for  the 
S^^\k^,  p,  5)  is,  therefore,  easily  obtained  by  analogy  to  the  solution  in  section  3.  In  calculating  the 
signal  from  this  solution,  we  use  one  more  approximation,  namely  we  replace  eq.  (4.14)  by 


s)  =  (2IN)  S  f>,  s) .  (4.23) 

We  thus  approximate  all  polariton  transformation  coefficients  for  the  exciton-like  polaritons  by 
unity.  Admittedly,  this  is  a  strong  discretization  of  the  exciton  character,  but  in  view  of  the  very  simple 
scattering  model,  it  seems  unnecessary  to  account  in  a  more  rigorous  way  for  these  coefficients.  The 
solution  for  W(/itg,s)  is  now  of  the  same  form  as  eq.  (3.46b),  with  y-*yp^  2nd 

^{p)^  %-k^i2  -  <^p*k^i2-  the  coherent  limit  (/j,  =  0),  this  yields  for  the  signal 

5(t)  =  exp(-27pT)  (coherent),  (4  24) 


which  agrees  with  our  general  conclusion  that  in  the  absence  of  scattering  the  signal  does  not  decay 
(apart  from  the  trivial  population  loss  with  rate  2yp).  In  the  opposite  limit  of  strong  scattering 
(fp  ^  Tp)  we  obtain,  in  analogy  to  eq.  (3.49), 

5(r)  =  exp[-2(yp  +  D^\k^  -  it2p)T|  (incoherent) ,  (4.25) 


with  the  polariton  diffusion  constant  (tc  isor) 
1  1 


en  ■; 


2  {(^p-jt^/2  ■ 


(4.26) 


Here,  we  assumed  that  ,  which  holds  under  general  symmetry  conditions  with  respect  to  the 

experimental  setup. 

Equation  (4.26)  can  be  rewritten  using  the  polariton  group  velocity  Up(p)  =  ^w^  [cf.  eq.  (3.38a)]. 
For  an  isotropic  rf-dimensional  system,  is  a  function  of  \p\  only  and  we  have 


(4.27) 
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with  an  effective  polariton  group  velocity.  As  polaritons  have  a  photon  component,  it  is  clear  that  c,. 
is  larger  than  the  effective  exciton  group  velocity  (eq.  3.38b)  for  the  same  system.  Thus,  ignoring 
possible  differences  in  the  exciton  and  polariton  scattering  rates,  the  general  effect  of  polaritons  on  the 
TG  will  be  to  accelerate  the  decay.  An  estimate  for  is  calculated  in  appendix  D.  within  the  infinite 
effective  exciton-mass  approximation  (J^  =  0  then!)  on  a  simple  cubic  lattice.  We  find 

(Up/c)-  =  TT(flO)(a/\,y .  •  (4.28) 


with  a  the  lattice  constant  and  A,,  =  lire  10  the  vacuum  wavelength  corresponding  to  the  molecular 
transition. 

Although,  strictly  speaking,  our  crystal  model  is  oversimplified  to  fully  represent  aromatic  molecular 
crystals,  it  is  tempting  to  use  eq.  (4.28)  also  for  this  class  of  crystals  by  replacing  a'  by  the  volume  of  a 
unit  cell.  We  then  obtain  17^  =  2  x  10^  cm/s  and  Up==10^cm/s  for  the  lowest  singlet  a-transitions  in 
naphthalene  and  anthracene,  respectively.  (For /and  (2.  we  used  the  values  given  in  section  4.1;  for  a  . 
we  used  400  A’  for  naphthalene  and  600  A^  in  the  case  of  anthracene.)  This  result  for  anthracene  is  an 
order  of  magnitude  larger  than  the  estimate  made  by  Agranovich  et  al.  [51]  for  the  relevant  polariton 
group  velocity  (see  below).  In  view  of  the  simple  model  and  the  approximations  used,  it  is  not  surprising 
to  find  such  a  large  discrepancy. 

The  situation  in  anthracene  is  probably  much  better  described  by  the  second  scattering  model  that 
we  wish  to  discuss  here:  the  bottleneck  model.  This  model  is  inspired  by  the  explanation  by  Agranovich 
et  al.  [51.  52]  for  the  TG  experiments  on  anthracene  crystals  carried  out  by  Rose  et  al.  [2b].  Let  the 
experiment  be  such  that  the  two  pump  pulses  have  frequencies  just  in  the  exciton  band,  so  that  they 
excite  first-order  polaritons  with  high  wave  vector  and  strong  exciton  character  (fig.  8).  Accounting  tor 
a  positive  effective  exciton  mass,  it  is  reasonable  to  assume  that  at  low  temperature,  polariton-phonon 
scattering  will  cause  the  initially  created  polaritons  to  relax  rapidly  to  the  bottleneck  region,  where  the 


Fig.  S  Illustration  of  the  bottleneck  irnxlel  for  ptilariton  scattering.  Compared  to  fig.  7,  a  larger  part  of  the  Bnllouin  zone  is  shown  ana  .1  finite 
(positive)  effective  exciton  mass  has  been  included.  High  wave-vector  polaritons  in  the  lower  branch  are  excited  by  pulses  with  frequencies  mst 
above  SI,  (the  zero  wave-vector  exciton  frequency)  and  rapidly  relax  to  the  polariton  bottleneck,  where  the  actual  TG  decay  takes  place. 
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density  of  states  levels  off.  This  classical  (particlef  picture  translates  into  the  following  model  for  rapid 
initial  relaxation  of  the  polariton  coherences  p,  t). 

(i)  The  polariton  center-of-mass  wave  vector  k  is  conserved  at  its  initial  value  =  -  k.:  we  have 

obtained  this  previously  as  a  general  result  for  phonon-induced  scattering  processes.  This  is  equivalent 
to  the  assumption  in  ret.  [52]  that  during  the  rapid  relaxation  process,  the  spatial  density  of  polaritons 
created  by  the  interfering  pump  pulses  is  unchanged. 

(ii)  p  relaxes  rapidly  to  the  bottleneck  region.  This  is  equivalent  to  the  statement  in  ref.  [52]  that  the 
polariton  momentum  (wave  vector)  relaxes  rapidly.  In  general,  p  must  be  expected  to  be  distributed 
over  the  entire  bottleneck  region  after  the  relaxation,  because  many  scattering  steps  involving  many 
different  phonons  can  take  place  during  this  process  [52].  We  note  that  the  classical  interpretation  of  p 
as  the  polariton  momentum  that  relaxes  to  the  bottleneck,  only  makes  sense  if  \kJ2\  is  small  compared 
to  the  wave-vector  width  of  the  bottleneck  region.  Only  then  does  the  coherence  p,  t)  involve 
two  polaritons  with  wave  vectors  p±  kJ2  that  also  lie  within  the  bottleneck  region.  For  the  experiments 
on  anthracene,  this  is  indeed  the  case,  as  these  were  carried  out  at  very  small  cross  angles  0  of  the 
pump  pulses,  so  that  the  grating  wave  number  was  always  much  less  than  optical  wave  numbers  [2b]. 

If  the  initial  relaxation  is  fast  enough,  the  observed  part  of  the  TG  decay  now  starts  from  a  new 
initial  condition  for  the  coherences  within  the  bottleneck.  We  will  assume  that  within  the  bottleneck  all 
polaritons  scatter  into  each  other  with  rate  /j,  and  that  there  is  an  overall  loss  of  polaritons  with  rate  y^. 
Like  the  first  scattering  model  that  we  discussed,  this  defines  a  restricted  Haken-Strobl  model,  for 
which  the  equation  of  motion  (4.16)  can  be  solved  analytically.  To  calculate  the  signal  from  this 
solution,  we  use  eq.  (4.23),  which  is  a  good  approximation  in  the  bottleneck,  because  the  coefficients 
jjCjJ  are  close  to  unity  there.  Nevertheless,  the  polariton  group  velocity  may  differ  appreciably  from  the 
bare  exciton  group  velocity  [60].  Let  denote  the  typical  group  velocity  in  the  bottleneck  region.  If 
Tp  >  IkglUh*  fhe  polariton  motion  is  diffusive  on  the  length  scale  of  the  grating  and  we  recover  eq.  (4.25) 
for  the  signal,  with  the  diffusion  constant  eq.  (4.27),  except  that  the  effective  sroup  velocity  i7p  is 
replaced  by  v^. 

This  result  is  independent  of  the  exact  initial  condition  within  the  bottleneck;  such  an  insensitivity  to 
the  initial  condition  was  also  noted  in  the  strong  scattering  limit  of  the  exciton  theory  (section  3.4).  The 
present  expression  for  the  signal  was  also  used  by  Agranovich  et  al.  [51.  52];  our  approach  shows  how 
to  derive  it  microscopically.  For  anthracene,  10^  cm/s;  combining  this  with  a  scattering  length 
vjtp  =  10”'  cm  (estimated  from  absorption  data),  one  arrives  at  Z>p  ~  1  cm7s  [51],  which  agrees  with 
the  observed  values  [2b]. 

We  conclude  the  discussion  of  the  bottleneck  model  by  noting  that  in  the  limit  of  weak  scattering. 
Tp  <  the  initial  condition  is  in  general  important  for  the  form  of  the  TG  signal.  Decay  will  always 
occur,  even  if  r^  =  y^=  0,  because  the  initially  populated  coherences  E{k^.  p.  /)  for  different  p  values  in 
the  bottleneck  will  dephase  with  respect  to  each  other. 

4.4.  Polariton  effects  in  frequency -domain  four-wave  mixing 

In  this  section  we  discuss  the  degenerate  four-wave  mixing  experiment  in  terms  of  polaritons.  Two 
cw  external  electric  fields  with  wave  vectors  it',  frequencies  a)-,  and  amplitudes  E]  (j=\.2)  are 
incident  on  the  crystal  and  we  are  interested  in  the  signal  with  wave  vector  k[  =  2k\  -  kl  and  frequency 
Wj  =  2a>,  -  o)^.  The  amplitude  of  this  signal  is  proportional  to  the  amplitude  of  the  polariton  with  wave 
vector  k^  =  2k^  -  k^  that  is  generated  nonlinearly  in  the  crystal  by  the  first-order  polaritons  with  wave 
vectors  k^  and  k^.  Our  description  is  based  on  the  temporal  Fourier  transforms  of  the  two-particle 
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polariton  hierarchy  eqs.  (4.13a)  and  (4.16).  We  wiir confine  ourselves  to  the  restricted  Haken-Strobl 
model  introduced  in  section  4.3,  and,  as  in  the  exciton  theory  of  section  3.3,  we  will  only  calculate 
contributions  to  the  signal  with  possible  resonances  at  oj,,  =  a>,  -  oj,  =  0. 

From  the  Fourier  transform  of  eq.  (4.13a)  we  find  for  the  signal  amplitude  to  third  order  in  the  field 
amplitudes 


^ (-f*  *'i)^t'!,|(<«^i) S p. w,,) , 


(4.29) 


where  we  used  eq.  (4.23)  and  k^  =  -  A,,  as  before.  The  two-particle  variables  for  different  momenta 

p  obey  coupled  equations  of  motion  that  are  obtained  from  eq.  (4.16)  with  the  restricted  Haken-Strobl 
scattering  model. 


-iw,.5'-’(*^,  p,  w,-)  =  [i( 


,  :)  -  +  rp)l-=‘'*(*g.  />.  ‘«>,:)  +  S  =''\k^,  p'.  w,,) 


+  t'p  *;)f ; + ^i;;(a>,  )A*(k,v,,  k',)Ei‘]8^ 


(4.30) 


p^^(k,+k,)/2. 


This  equation  is  solved  in  the  standard  way  (appendix  B)  and  after  substituting  the  first-order  polariton 
amplitudes  from  the  linear  approximation  to  eq.  (4.13a),  we  obtain 

.  y{k,u,)A%u,A[)A*(k.v,,k:)iE'l)-E^.* 

^  : - - - - — t; - 

-ita,,  -  i(ca*,  -  w*  )  +  Tp  +  Yp 


/  V  •  • 

V~  ^  (-I"i2  -1%-*. ;  +  ^  r  +  y  y 


We  will  assume  that  the  second  factor  on  the  r.h.s.  has  no  sharp  dependence  on  w,-,.  For  excitons.  this 
was  found  as  a  result  of  cancellations  between  the  source  cr,^  and  the  quantity  -ica,,  -  i[y(itO  ~ 
•^(^i)l  +  /’  +  y  in  eq.  (3.32).  In  the  present  case,  io-^  is  replaced  by  the  factor  multiplying  5^^  in  eq. 

(4.30).  and  a  cancellation  is  not  likely,  even  though  we  have  no  detailed  form  for  A(kv.k')-  Note, 
however,  that  w,  =  and  <a,  =  ,  as  the  frequencies  of  the  first-order  polaritons  have  to  equal  the 

frequencies  of  the  external  fields.  This  is  an  essential  difference  with  the  exciton  theorv. 

We  thus  conclude  that  the  first  denominator  on  the  r.h.s.  of  eq.  (4.31)  does  not  depend  on  w, ,  at  all! 
Any  w,, -resonance  must  emerge  from  the  last  factor. 


R{k^,  w,.)  =  (l  -  E 


^m-ir  *» 

P  - 


+ ^P  +  Tp)"')  • 


(4.32) 


For  Tp  0,  we  have  R{k^,  ta,,)  -  1.  so  that  no  ta,, -resonance  is  found.  Conversely,  in  the  limit  of  strong 
dephasing  y^,  ~  2I)’  (cf-  appendix  C), 

R{k^,  (a,,)  =  1  +  rp/(-i(a,,  +  Yp  +  Dp*;) . 


(4.33) 
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with  £)p  the  polariton  diffusion  constant  as  defined  in  eq.  (4.26).  The  intensity  of  the  D4WM  signal, 
a>p)|',  then  exhibits  a  Lorentzian  resonance  analogous  to  the  one  in  eq.  (3.35) 
(replace  t-*  F^,  y— >  yp,  Dp).  Comparing  this  with  our  results  in  section  4.3,  we  also  conclude 

that,  as  in  the  exciton  case,  the  amplitude  of  the  D4WM  signal  is  related  to  the  amplitude  of  the  TG 
signal  by  a  single  Fourier  transform. 

5.  Concluding  remarks 

The  calculation  of  the  nonlinear  optical  response  in  condensed  phases  is  a  complicated  many-body 
problem  involving  the  material  and  the  radiation-field  degrees  of  freedom.  Numerous  approximate 
schemes  have  been  developed  over  the  years  with  various  degrees  of  sophistication  in  order  to  calculate 
optical  nonlinearities  [19-23,  70,  71,  92].  The  present  review  provides  a  unified  theoretical  description 
which  clarifies  the  interrelationships  among  the  different  schemes  and  their  range  of  validity.  Figure  9 
illustrates  the  systematic  hierarchy  of  approximations.  Starting  with  the  multipolar  Hamiltonian  [73] 
which  is  fully  retarded  and  contains  no  explicit  intermolecular  interactions,  we  can  rigorously  derive 


Fig.  9.  Schematic  diagram  showing  the  hierarchy  of  approximations  used  in  this  paper. 
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Heisenberg  equations  of  motion  (eq.  2.9)  in  which  Tntermolecular  forces  are  explicitly  included  [76|. 
Equation  (2.9).  which  applies  to  any  material  and  Held  operator  Q,  is  the  starting  point  for  the  present 
microscopic  theory.  For  Q  =  E,  we  get  the  Maxwell  equation  (cq.  4.3b).  In  simple  single-molecule 
theories  of  nonlinear  response,  one  often  uses  the  Schrodinger  picture  and  solves  for  the  density  matrix 
(19-23).  This  procedure  is  impractical  for  complex  many-body  systems,  since  the  complete  many-body 
density  matrix  is  very  complicated  and  contains  too  much  information.  The  Heisenberg  equations,  on 
the  other  hand,  allow  the  direct  calculation  of  the  relevant  observables. 

At  this  point,  as  shown  in  fig.  9,  we  should  follow  two  separate  routes,  depending  on  the 
radiation-matter  coupling  strength  (/).  The  key  parameter  is  /7r(jk).  The  exciton  dephasmg  rate  r(k) 
may  be  obtained  from  linear  optical  measurements  [see  eq.  (3.12b)].  If  it  is  sufficiently  large ]r(fc)  ^/j. 
it  destroys  the  coherence  between  the  radiation  and  e.xciton  modes,  which  then  become  weaklv 
coupled.  In  this  case,  the  material-radiation  correlations  can  be  neglected,  and  we  mav  derive  coupled 
equations  for  the  field  and  the  material  degrees  of  freedom  as  illustrated  in  the  left  route  in  fig.  9 
(section  3).  The  field  satisfies  the  Maxwell  equations  whereas  the  material  system  (excitons)  satisfies  an 
infinite  hierarchy  of  coupled  equations  in  which  the  polarization  (/*)  is  coupled  to  the  expectation  value 
of  the  Maxwell  field  and  to  higher  variables  such  as  (PW),  which  in  turn  couple  successively  to 
higher-order  variables.  By  solving  these  equations  perturbativelv  in  the  average  Ma.xwell  field  E.  we 
obtain  the  optical  susceptibilities  A'*"’  t'  '-  <^tc. 

The  success  of  this  method  is  based  on  the  fact  that  the  linear  optical  properties  of  the  system  depend 
only  on  single-particle  states  whereas  weak  nonlinearities,  such  as  and  depend  only  on  a  few 
particles.  This  situation,  which  is  similar  to  the  zero-temperature  many-body  theory  [44.  81).  allows  us 
to  truncate  the  hierarchy  very  early  and  still  maintain  the  essential  physics  of  the  system.  The  simplest 
approximation  is  obtained  by  factorizing  at  the  single-particle  level  (section  3.2).  In  this  case,  the  only 
relevant  material  dynamical  variables  (by  assumption)  are  the  polarization  variables.  The  ■'onlinear 
susceptibility  y'''  can  then  be  written  ih  various  equivalent  forms.  Within  this  factorization,  the 
local-field  approximation  (eqs.  (3.19)  and  (3,20)[  holds,  and  the  many-body  problem  reduces  to 
essentially  a  single-body  calculation.  We  further  recover  the  coupled-oscillator  expression  (eq.  3.14a) 
with  interaction-induced  correction  terms.  At  this  level  of  description,  the  susceptibilities  have 
resonances  at  the  Frenkel  exciton  energies,  cxciton  transport  is  not  accounted  for.  and  we  cannot 
describe  the  TG  or  degenerate  four-wave  mixing  experiments. 

When  the  hierarchy  is  truncated  at  the  binary  (two-particle)  level  (sections  3.3  and  3.4).  we  obtain  a 
good  description  of  quantum  transport  and  we  can  make  the  connection  to  the  phenomenological 
macroscopic  treatments  of  the  grating  experiment  which  use  master  equations,  the  Boltzmann  equation, 
and  the  diffusion  equation  (section  3.5).  The  single  Fourier-transform  relation  between  the  TG  and  its 
frequency-domain  analogue  (degenerate  four-wave  mixing)  is  demonstrated. 

The  two-particle  description  illustrates  the  limitations  of  Bloembcrgen  s  classical  anharmonic  oscil¬ 
lator  picture  [19].  When  nonlinearities  are  important,  new  dynamical  variables  (other  than  the 
polarization  P)  become  relevant.  The  oscillator  model  neglects  the  contributions  of  these  additional 
variables  and  requires  the  factorization  of  every  variable  in  terms  of  the  polarization.  In  our  equations 
we  choose  a  particular  two-body  factorization  in  which  \V  is  the  source  of  nonlinearity.  This 
nonlinearity  is  local  in  space  (intramolecular)  and  seems  adequate  for  the  calculation  of  the  TG 
experiment,  since  the  signal  depends  on  VV.  It  should  be  emphasized,  however,  that  alternative  sources 
of  nonlinear  evolution  may  be  important  in  other  situations  [47,  59]  (cf.  section  3.1).  In  particular, 
enhanced  (cooperative)  optical  nonlinearities,  which  may  be  observed  in  molecular  aggregates,  have 
been  shown  to  originate  from  intermolecular  (nonlocal)  nonlinearities  [35].  The  present  hierarchy  is 
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illustrated  in  fig.  9.  The  single-particle  description  can  be  obtained  directly  from  eqs.  (2.9)  (as  done  in 
section  3.2).  or  by  factorizing  the  two-particle  variables  into  products  of  single-particle  variables. 

When  the  oscillator  strength  per  unit  volume  is  sufficiently  large  and  dephasing  processes,  are  slow 
[/^  the  radiation  and  the  exciton  degrees  of  freedom  are  strongly  coupled,  and  they  should  be 
treated  as  a  single  dynamical  system.  In  this  limit,  the  description  of  optical  processes  in  terms  of 
polariton  models  is  most  natural  (sectiorr  4)  and  is  illustrated  in  the  route  shown  on  the  right  side  in  fig. 

-  9.  Although  it  is  formally  possible  to  define  nonlinear  susceptibilities  even  in  this  limit  (66.  67].  they  are 
complicated  and  not  particularly  useful.  In  section  4  we  develop  a  polariton  hierarchy  [76]  which  allows 
us  to  expand  the  nonlinear  signal  directly  in  terms  of  the  external  fields  (rather  than  the  Maxwell  fields) 
and  to  abandon  the  notion  of  susceptibilities.  Note  that  this  expansion  then  depends  on  the  boundarv 
conditions  relating  the  external  field  to  the  polaritons. 

Ovander  [59]  has  incorporated  polariton  effects  in  nonlinear  optics,  by  expressing  his  total  Hamil¬ 
tonian.  which  also  contains  cubic  and  quartic  terms,  in  terms  of  the  polariton  creation  and  annihilation 
operators.  The  harmonic  (quadratic)  part  then  assumes,  by  definition,  the  simple  form  of  a  set  of 
noninteracting  harmonic  oscillators  (the  bare  polaritons).  The  anharmonic  terms  yield,  even  for  the 
simplest  approximations  of  the  polariton  transformation,  a  large  number  of  interaction  terms  which  are 
responsible  for  nonlinear  processes,  and  Ovander  treats  them  as  perturbations  to  the  harmonic 
Hamiltonian.  The  Fermi  golden  rule  is  used  to  evaluate  polariton  fusion  and  fission  rates,  which  are 
related  to,  e.g.,  the  sum  harmonic  intensity,  the  Raman  scattering  intensity,  etc.,  depending  on  the 
specific  perturbation  term  under  consideration.  This  formulation  has  several  limitations.  The  resulting 
Hamiltonian  is  very  complicated,  the  Bose  approximation  used  for  the  Pauli  operators  results  in  the 
neglect  of  intramolecular  nonlinearities,  ar  '  the  Fermi  golden  rule  does  not  suffice  to  treat  nonlinear 
processes  of  order  higher  than  three:  higher-order  perturbation  theory  has  to  be  used  to  describe  these. 

Our  polariton  hierarchy  does  not  suffer  from  these  difficulties.  In  the  strong-coupling  regime 
analyzed  in  section  4,  the  TG  probes  the  dynamics  of  polaritons.  Like  the  exciton  motion,  polariton 
dynamics  can  be  either  coherent  or  incoherent  (diffusive),  depending  on  the  magnitude  of  polariton 
dephasing.  Two  conditions  have  to  be  satisfied  in  order  for  the  dynamics  to  show  diffusive  polariton 
motion,  namely.  Fik) « /and  t  ^  |w-_*  ;  -  w  *  ,|  (section  4.3).  The  exciton  dephasing  rate  has  to  be 
sufficiently  small  for  the  elementary  excitations  to  be  polaritons  and  the  polariton  dephasing  should  be 
sufficiently,  large  to  make  the  motion  incoherent.  Both  conditions  can  be  satisfied  simultaneously.  A 
two-particle  and  single-particle  level  of  approximation  exists  for  polaritons.  in  complete  analogy  with 
what  we  derived  for  excitons.  Since  our  interest  is  in  polariton  transport,  which  is  absent  in  the 
single-particle  description,  we  considered  only  the  two-particle  level. 

An  important  conclusion  of  the  present  anaivsis  is  that  observation  of  coherent  exciton  motion  for 
dipolar  excitons  is  impossible,  because,  when  exciton  damping  is  small  enough  to  give  coherent  exciton 
motion,  the  radiation  modes  are  necessarily  strongly  coupled  with  the  polarization,  so  that  polariton 
effects  become  significant.  Coherent  exciton  motion  over  distances  greater  than  the  optical  wavelength 
is  therefore  an  unrealistic  and  oversimplified  model  of  elementary  excitations  in  molecular  crystals  at 
low  temperatures. 

The  present  formulation  and  equations  of  motion  can  be  directly  used  and  extended  to  treat  optical 
nonlinearities  in  other  systems.  There  is  currently  a  growing  interest  in  the  optical  properties  of 
nanostructures  [108-113].  These  are  fabricated  molecular  assemblies  with  specific  molecular-level 
order.  Examples  are  clusters,  monolayers  and  multilayers.  Our  equations  of  motion  are  ideally  suited 
for  the  microscopic  calculations  of  optical  nonlinearities  in  these  systems.  Other  interesting  extensions 
of  the  present  formulation  involve  the  incorporation  of  short-range  forces  and  static  and  dynamical 
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disorder  [82].  This  is  important  for  calculating  optical  nonlinearities  in  doped  glasses,  polymers, 
solutions,  and  in  mixed  crystals. 

Acknowledgements 

J.K.  acknowledges  stipport  by  the  Netherlands  Organization  for  Scientific  Research  (NWO)  through 
the  award  of  a  Huygens  Fellowship.  S.M.  acknowledges  the  Camille  and  Henry  Dreyfus  Foundation. 
The  support  of  the  Air  Force  Office  of  Scientific  Research,  the  National  Science  Foundation  and  the 
Petroleum  Research  Fund,  administered  by  the  American  Chemical  Society,  is  gratefully  acknowl¬ 
edged.  We  wish  to  thank  Professor  V.M.  Agranovich  for  useful  comments.  Finally,  we  thank  Rose 
Marie  Ferreri  and  Arlene  Bristol  for  the  careful  typing  of  the  manuscript. 


Appendix  A.  Elimination  of  phonon  degrees  of  freedom 

In  this  appendix  we  discuss  the  elimination  of  phonon  variables  underlying  the  various  self-energy 
and  scattering  kernels  for  excitons  and  polaritons  given  in  sections  3  and  4,  respectively.  To  this  end. 
several  approaches  can  be  used,  which  do  not  differ  essentially  in  their  results  [42.  bS.  114).  Here  we 
use  the  projection-operator  technique  which  was  also  employed  in  ref.  [68].  The  relevant  Hamiltonian 
for  our  problem  reads 

IA,n 

Here  stands  for  the  "system”  Hamiltonian.  To  describe  the  scattering  of  excitons  on  phonons 
(section  3).  we  choose  the  exciton  Hamiltonian  for  H^. 

H,  =  ft  2  .  (A.;al 

k 

(We  use  the  Heitler-London  approximation  in  this  appendix.)  To  describe  polariton-phonon  scattering 
(section  4).  we  choose 


/i  ^  • 


(A.:b) 


with  Wj,,  the  polariton  dispersion  in  the  ab.sencc  of  exciton  damping  [/’(A)  =  ()  in  eq.  (4.b)].  Both 
excitons  and  polaritons  are  treated  as  bosons  in  the  calculation  of  the  self-energies.  The  phonons 
constitute  a  "bath”  and  are  described  by  the  Hamiltonian  (^4-  3.13).  Finally.  denotes  the 
system-bath  interaction  (eq.  2.14).  If  we  take  the  polaritons  as  the  system.  has  to  be  expressed  in 
terms  of  polariton  operators,  which  is  done  using  the  inverse  polariton  transformation  eq.  (4.bb);  it  is 
customary  to  make  the  approximation  c*,,  =0  when  doing  this  ]60a.  68.  80]. 

We  will  work  in  the  Schrodinger  picture  and  define  p{t)  to  be  the  total  density  operator  of  the  system 
and  the  bath.  We  now  define  a  projection  onto  a  reduced  density  operator  by  [68.  1 1.''] 


48 


J  Knoesier  ami  S  Mukamel.  Transient  gratings,  four-nave  mixing  and  polariion  effects  in  nonlinear  optics 


where  is  the  density  operator  for  phonons  in  tTtermal  equilibrium  at  temperature  T  and  Tr^„  traces 
over  the  phonon  degrees  of  freedom.  The  details  of  this  projection  can  be  found  in  ref.  [nS|.  In 
particular,  it  can  be  shown  that  up  to  second  order  in  the  contiibuiion  from  the  phonon  bath  to  the 
evolution  of  the  expectation  value  of  an  arbitrary  system  operator  O^.  is  given  by 

l 

(  ^  =  -  J  d/'Tr[{L,pexp[ii4i-/')]  o(/')] .  (A. 4) 

—  X. 

=  =  (A.Ja.b) 


where  Tr  denotes  the  total  trace. 

We  first  consider  the  exciton  case  and -calculate  the  effect  of  the  phonon  bath  on  the  variables 
(B^fO)  and  .  For  {B^{t)),  the  result  is  readily  obtained  from  the  scattering  contribution 

that  was  derived  in  appendix  A  of  ref.  [68].  To  this  end.  we  drop  the  branch-labels  and 
-sums  in  that  derivation,  we  replace  the  polariton  dispersion  by  the  e.xciton  dispersion  and  we  set  all 
polariton  transformation  coefficients  equal  to  unity.  Without  further  derivation  we  give  the  result 
within  the  Markov  approximation.,  keeping  both  the  frequency  shift  and  width. 

[(d/dO(B*(/))U  =  -iIU)(B*(t)),  (A.6a) 

with  the  complex  self-energy 

_ Ot  \ 

(A.bb) 

f\(k.  <7)  =  q)  +  x,(q)  . 

Here  ) ,  is  the  thermal  (Bose- Einstein)  occupation  of  phonons  with  wave  number  q  in  branch  5. 
We  note  that  this  result  can  also  be  derived  on  the  operator  level  [114],  which  we  used  in  cq.  (3.1). 
Taking  the  limit  77— »()  explicitly,  we  arrive  at  eqs.  (3.3)  for  the  phonon-induced  frequency  shift  and 
damping  rate  of  the  exciton.  Of  course,  {B^{i))  obeys  eq.  (A. 6a)  with  2.(jt)  replaced  bv  -l*(jt) 
[2  *(A:)  denotes  the  complex  conjugate  of  2(it)j. 

We  now  bnefly  outline  the  derivation  of  the  phonon  contribution  to  the  evolution  of  the  exciton 
coherence  {B^{t)B^(t)),  which  we  need  in  section  3.3.  We  take  d^  =  B^B^  in  eq.  (A. 4).  and 
straightforward  algebra  leads  to  (t  =  /  -  t') 


'it;;  2  \FSk,  ^)|- 


h  N  r|— *0' 


^k+q  fi  irj 


T  pexp(iL,T)  L^^(B,B, 


^  I  9)exp[i(/2,.^  -  Fi,  )rl  J',^(7)| 

-  E  F^{k’  -  q,  q)c\p[\{n^  -  Fi*  _,)t]  L^^[BlB,  _^X^^(r]\  . 

(A. 7a) 


exp(-ir},,-)  +  exp(i/7^,,T) . 


(A. 7b) 
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For  the  commutators  appearing  in  the  r.h.s.  of  eq.lA.7a)  we  obtain 

"hhj:  1^'*  ~  ^  ■  «'■  1'  's*  .  H*. .  - 

1 
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(A. 8a) 


^-,  v  ) 

(A.8b) 

We  will  neglect  the  last  contributions  (the  k"  summations)  in  these  equations,  as  ihcv  couple  to  a  higher 
(quartic)  exciton  variable.  Wc  now  combine  cqs.  (.A. 7)  and  (A. 8).  and  after  stime  straightforward 
algebra  we  obtain 

Tr({Z-,pexp(iL^7)  L^^{BlB^  )]'&(t')\=^  -i-  ^  f  ;(fc.  q)exp|-i(fi*  -  n^.)7]G{T) 

n  ;V  V' 

^{n(k  +  q--q){K(nB,(t’)}-F:(k'  +  q.-q){B,  .y)B,.^(n)] 

-  ^  Z  P(*'  +  q.  -</)exp(-i(fi,  -  (}^  ,^)r\G^(r) 

X  [F:{k.  q){K-,(nBx^y))  -  p(a-'.  q){B,  {nB,(i'))\.  (a.9) 

(«,>)rexp(if\,T)  +  +  l),.exp(-if}  ^,7).  (A. 10) 

Equation  (A.9)  combined  with  oq.  (A. 4)  defines  the  phonon  contribution  to  the  equation  of  motion 
tor  (/)fij(0).  This  equation  can  be  made  time-local  by  applying  a  Markov  approximation.  To  this 
end.  we  approximate  [68 1 

(fi;  (/')fi,(/'))  ^{Bl  (t)B,{t))  exp[-i(/4  -  /2j7| .  (A.ll) 

and  we  calculate  the  remaining  time  integrals  by  adding  a  convergence  factor  expf-T^r)  (77— »0  ).  We 
then  obtain 

[(d/do(fi*  =  i{iHk')-  :iik)\{B,  (i)B,(i)} 

+  i  Z \^*ik\ k.  q)  -  i(*.  k'.  9)1(  s,.  .,(r)B,.,(r))  . 


-  A(t)] 

=  -j^.  Z  {F,(k.  q’)B[^^  fl,  -  F,(k'  -q-q'.  q')B,B,.  , 

+  Z  F^{k'\  -Jexp(i/7  ,,7)-exp(-i/7,^r)] . 


(A.i:) 
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with  I{k)  as  defined  in  eq.  (A. 6b)  and 


(^v)T 


o. 


+ 


(«-,v  +  l)l 


~  1^7  -  117 


(A.13) 


The  first  r.h.s.  term  in  eq.  (A. 12)  simply  contains  the  uncorrelated  sum  of  the  self-energies  of  {B^  (t)) 
and  and  is  the  “ 7, "-contribution  [22.  75)  to  the  total  self-energy  of  {Bl.{t)Bi^{t)) .  The  second 

r.h.s.  term  in  eq.  (A.  12)  reflects  correlated  dynamics  of  the  bra  and  the  ket  side  of  the  density  operator 
(“7*”  or  pure  dephasing  contribution)  [22.  75).  If  we  rewrite  eq.  (A.  12)  in  terms  of  the  variables 
Q(k,  p,  t)  introduced  in  section  3.3,  we  finally  arrive  at  the  phonon  contributions  in  eq.  (3.25)  with  the 
self-energy  matrix  eq.  (3.26). 

We  conclude  this  appendix  by  considering  the  case  that  the  polaritons  constitute  the  system  (eq. 
A. 2b)  and  we  give  the  phonon-induced  self-energies  for  (i*^(7))  and  In  the  approxi¬ 

mation  that  the  polariton  transformation  coefficients  =  0,  the  calculation  is  very  analogous  to  the 
one  for  excitons.  and  we  just  give  the  final  results. 


\d 

.dr 


(L(0) 


x.Ak)a,At)). 


(A.  14a) 


^vAk)^-  lim  S  2 \F,{k. 

fl  1\  T7-»0  qs 


-I- 


+  I)l 


(X), 


k  +  qv" 


~  %  (^k-q.  -  <^kv  +  '  >^7 


(A.  14b) 


lj;  {(iAOLU)) 


ph 


u  r  " 

+  '22  (i*.. k. ,)  -  q)i((:-,.AnL,.(0>  ■ 


,.4*.  *■.  ,)=  ^  Urn  2  n(k.  q)F,(k\  q)x„x:.,,  x:.,  x,  -...- 


( 


("„)t 


{n  qs  +  Or 


^kxqv 


^k  v  -  iiqs  -  f^k  -q.  -  -  W*  ,,  +  ii-qx  "  'V 


(A. 15a) 


(A.  15b) 


Note  that  Eqs.  (A.  14b)  and  (A.  15b)  reduce  to  the  e.xciton  self-energies  eqs.  (A. 6b)  and  (A.13). 
respectively,  if  all  are  set  to  unity,  the  (summations  over)  branch  labels  are  dropped,  and  the 
polariton  dispersion  is  replaced  everywhere  by  the  exciton  dispersion.  Equation  (A.  15a)  for  the 
two-particle  polariton  variables  can  be  rewritten  in  terms  of  the  variables  E,.„  ik.  p.  i)  introduced  in 
section  4.2.  This  leads  to  the  second  term  in  eq.  (4.16).  with 


Ak.  ix:  p,p')^  [8^.„I,Ap  +  k/2)  -  k/2)]8^^  +  I,„  „  „  (p  +  k/2.  p  -  k2.  p  -  p) 


_  V 


,-,,-Ap-  k/2.  p  +  k/2,  p'  -p). 


(A.16) 


where  v  is  symbolic  for  (r.  v'.  v".  u'"). 
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Appendix  B.  The  two-particle  Green  function 

In  this  appendix,  we  solve  the  equation  of  motion  for  two-particle  coherences  within  the  Haken- 
Strobl  model.  This  equation  is  encountered  in  sections  3.3  (eq.  3.30).  3.4  (eq.  3.43),  4.3  (eq.  4.22)  and 
4.4  (eq.  4.30).  For  explicitness,  we  adhere  to  the  problem  of  section  3.3;  the  other  applications  are 

readily  translated  to  the  notation  used  here-.-  The  equation  of  interest  reads 

-  -  -  - 

p,  !«)  =  |i.?(  p)  -  (r + y  )|Q(/,,  «,)  +  (f/N)  1 0(  p'.  <«) + /( P) .  I B.  I ) 

p 

where  we  defined 

.fip)^J{p-  kJ2)  -  y(p  +  kjl) .  (B.2) 

The  variable  k^  is  suppressed  in  all  quantities,  as  the  equation  of  motion  is  trivial  with  respect  to  it  (A\.  is 
conserved).  The  oj-dependence  is  indicated  explicitly,  as  this  variable  is  independent  of  k^  in  the 
time-resolved  applications  (sections  3.4  and  4.3).  Finally,  we  assumed  a  general  source  f{p). 

For  one-dimensional  systems,  eq.  (B.l)  has  been  solved  by  Haken  and  Strobl  [88],  Reineker  [54], 
and  Garrity  and  Skinner  [55],  using  eigenvector  analysis.  In  this  appendix,  we  will  solve  eq.  (B.l)  using 
the  Green  function  method  of  ref.  [40].  The  problem  in  the  equation  is  the  second  right-hand  side  term, 
which  results  from  pure  dephasing  processes  (~0.  In  the  absence  of  this  term,  the  equation  is  diagonal 
in  p  and  trivial  to  solve: 

Q\p.  (.) -  G';(a;)/(p) .  Gl(to)  =  [-i<0  -i.^(p)+r+yr'.  (B.3a. b) 

and  expression  (B.3b)  is  the  unperturbed  Green  function. 

As  the  perturbation  in  the  full  equation  (B.l)  couples  to  all  momenta  p'  with  equal  strength,  it 
represents  a  single  impurity  in  the  Fourier  space  conjugate  to  p.  Therefore,  the  full  problem  can  be 
solved  through  a  T-matrix  analysis  [116].  We  use  the  following  Fourier  transforms: 

Q{m,  w)  =  ( 1  IVN)  2  Q{  p-  w)  exp(-ip  •  r^) ,  ( B.4a) 

p 

=  ^'ZHp)  exp[  -  ip  •  ( r,„  -  r  J] .  /(m)  =  f{  p)  exp(  -  ip  •  r,„ ).  ( B  .4b.  c) 

Equation  (B.l)  then  translates  into 

2  -  '<* . +(/■+  -  fS„.S„ole(".  =  <B.5) 

n 

The  full  Green  function  in  coordinate  representation  is  defined  by 

Q{m.to)^lG„Jto)f{n),  (B.6) 

n 

and  is  now  easily  found  to  obey 

S  [-iw5,„„,.  -  ii,„„,  +  (r  +  y)5„„.  -  r8„„  S„„]G„  „{(o)  =  .  (B.7) 

m ' 
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Standard  r-matrix  analysis  with  as  perturbing  potential  now  yields  [40] 

+  Cl(<»)G;L(«.)n|l  -  /•c:,(<.)| .  (B.8) 

Here,  is  the  unperturbed  Green  function  in  coordinate  representation 

Cl(«')  =  (l/N)Sc'>)exp|-i,,-(r„-r„)|  (B.9) 

P 

[cf.  eq.  (B.3b)].  The  eventual  quantity  of  interest  is  (see,  e.g..  eq.  (3.29)]  Z^Qip.  oj).  which  can  be 
written  as: 

iQip,  a;)  =  VNQ{m=Q,  to)  =  Vn{1  G:Ua,)/(/i))[l  -  tGlMV'  .  (B.lO) 

p  'n 

where  in  the  last  step  we  used  eqs.  (B.6)  and  (B.8).  Using  eqs.  (B.3b),  (B.4c),  and  (B.9),  we  finally 
arrive  at 


Yj  Q{p,  w)  = 


L^[-m-iJip)  +  r+y\-'f{p) 

i-(f//v)E,,[-ica-ii(p')  +  r+rr'  ■ 


(B.ll) 


Appendix  C.  The  diffusive  limit  of  the  D4WM  intensity 


In  this  appendix  we  derive  eq.  (3.35)  for  the  D4WM  intensity  in  the  diffusive  limit  [r>\J{p  -  k^l2) 
-J{p  +  k^l2)\  (for  all  p)  and  r>  y\.  We  start  from 


S(k^,  (oJ--^\R{k^,  ca,2)|‘  , 

(C.l) 

/?(*g,a>p)  =  (l-^  2[-it«^,2-i>^(p)+/’  +rr')  ’ 

(C.2a) 

where  we  used  the  shorthand  notation 

.9{p)  =  J{p-kJ2)-J{p  +  kJ2). 

(C.2b) 

Straightforward  Taylor  expansion  yields 

Y  [-i‘«>i:  -  '^HP)  +  r+y\~' 


-iwp  +  r+  y 


V  1 

Y  1  +  ~ — 

p  G  —  IWp 


i-^(p) 
p  +  r+  y 


•Hp) 

-ia>p  +  T  +  y 


) 

(C.3) 


The  second  term  on  the  r.h.s.  vanishes  upon  summation  [T.pJ{  p)  =  0],  because  J'{p)  is  the  difference  of 
two  lattice  Fourier  transforms  and  p  runs  over  an  entire  Brillouin  zone.  In  the  third  term  within  the 
r.h.s.  summation,  we  approximate  -iwp  +  r+ y  =  T,  assuming  that  over  the  frequency  range  of 
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interest  This  will  be  verified  a  posteriori.  Substituting  eq.  (C.3)  into  eq.  (C.2a).  we  then 

arrive  at 


r./t  .  ~iwp  +  F  +  y  F 

R{k  ,  iOp)  1=1  + 

-iwp  +  y  +  -iwp  +  y  + 

(C.4) 

(C.5) 

P 


At  this  point,  it  is  useful  to  note  that  R{ks  tu)  is  essentially  the  Green  function  for  the  exciton 
population  W{k,  to)  in  momentum-frequency  representation.  This  is  seen  by  combining  eqs.  (3.24b). 
(3.32),  and  (C.2a).  In  the  last  form  of  eq.  (C.4),  the  term  1  may  be  neglected  relative  to  the  second 
term,  so  that  we  recognize  the  Green  function  for  a  diffusing  particle  with  diffusion  constant  D^..  Thus. 
Dj  has  the  meaning  of  an  exciton  diffusion  coefficient;  in  general.  D,  depends  on  the  wave  vector  A:,,. 
Combining  eqs.  (C.l)  and  (C.4)  and  using  t  (because  |3f(p)|  F  for  all  p),  we  finally  obtain 

S{k,,  coJ  =  I  +  Fit  +  2y) /[to],  ^iy  +  D,k;f] .  (C.6) 

This  signal  has  a  Lorentzian  resonance  with  width  y  +  D^k'^  F,  which  justifies  a  posteriori  our 
assumption  that  Iw,,]  T  in  the  frequency  region  of  interest. 


Appendix  D.  Effective  polariton  group  velocity  in  the  restricted  Haken-StrobI  model 

In  this  appendix,  we  calculate  an  estimate  to  the  effective  polariton  group  velocity  defined  through 
eq.  (4.27)  for  the  restricted  Haken-Strobl  model.  Our  system  is  not  isotropic  and  the  polariton 
frequency  depends  both  on  the  direction  and  the  magnitude  of  its  wave  vector.  In  order  to  obtain  a 
simple  analytical  result,  we  will  replace  the  dispersion  relation  eq.  (4.6a)  [with  F{k)  =  ():  strong 
polariton  effects]  by  that  of  transverse  polaritons  in  an  isotropic  (atomic)  system  with  transition  dipole 
fi.  This  dispersion  relation  reads  {to  is  the  polariton  frequency)  [43.  48a.  57] 

{kcyid)'  =  1  +  /'/(at"  -  w")  =  tfw) ,  (D.l) 

with /■  =  HTrpiifi'lh  [cf.  eq.  (4.6b)].  As  mentioned  in  section  4.1./  is  the  frequency  separation  between 
the  upper  and  lower  polariton  curves  at  the  wave  number  A:,,  where  the  exciton  band  and  photon 
dispersion  line  cross.  In  practice, /2.  In  eq.  (D.l)  an  infinite  effective  exciton  mass  is  assumed  and 
(t)^  is  the  transverse  dipolar  exciton  frequency  at  optical  wave  vectors  (k~0)  ]48a]. 

-f-13)'  -^n-f-ibn .  (D.2) 

We  also  define  the  longitudinal  exciton  frequency  (48a), 

w,,  =  {fl-  +  2/73)''^  =  +  /'/3/2  .  (D.3) 

These  expressions  for  to^  and  W||  are  easily  obtained  from  eq.  (3.6)  and  the  approximation  to  the 
dipole  sum  J(k)  given  below  (eq.  3.17).  In  the  continuum  limit  for  k  (A.  N/V=  p).  we  now  have 
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(D.4) 


The  integral  has  two  contributions:  one  from  the  upper  branch,  extending  from  ^  =  0  to  /t  =  k„,  and  one 
from  the  lower,  ranging  from  k  =  k^^  to  the  edge  of  the  Brillouin  zone  (fig.  7).  It  is  convenient  to 
transform  to  the  frequency  as  variable. 


dt  =  j 


do)  (oi) . 


(D.5) 


Up(a))  can  be  found  by  taking  differentials  of  both  sides  of  the  relation 
e((v)a>'  =  k'C  . 


(D.6) 


We  then  obtain 


(D.7) 


Using  eq.  (D.l).  we  arrive  at 


Up(a))  =  cVf(w) 


{u)'^  -  10')'  +  /■&»■ 


(D.8) 


The  V^(<v)  poses  no  problem  as  the  integration  occurs  outside  the  region  where  5(a>)  <  0  (the  polariton 
stopgap  ranging  from  (o^  to  Wn). 

Combining  eqs.  (D.1)-(D.5)  and  (D.8).  we  arrive  at 


=  +  /, . 


u.,,  f : 


=  / 
477  pc  ■> 


d(U 


I  T  T  -* 

(0'((0'  -  Wjj)  '{(O'  -  CO'  ) 

T  T  •»  •)  1 

{(O'  -  (O'  )■  +  f'(0' 


:  j  : 


(D.9a) 


(D.9b) 


lO 


(0'{(0n-(0  ‘(W*  -  (O') 


{(O'  —  (O'  )'  +  f'(0' 


(D.9c) 


Here  we  have  used  the  fact  that  the  upper  branch  has  as  lowest  frequency  (o.  and  the  lower  branch  has 
highest  frequency  a>  .  The  other  integration  boundaries  are  determined  by  the  polariton  frequencies  at 
k  =  (o  =  (o  ±//2.  where  in  eq.  (D.9b)  we  approximated  (o  +//2  =  w,,  +/(1  -  flQ)ll~  (On  +//2. 
(Use  eqs.  (D.2)  and  (D.3)  and  the  fact  that  I.)  As  the  integration  domains  are  restricted  to  a 
region  close  to  the  transition  frequency  /3.  a  resonance  approximation  can  be  applied  to  the  integrands. 
After  straightforward  changes  of  variables,  we  arrive  at 
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n  f 


/,  =  -^  I  dx 


x'  ■(  1  +  .v) 


1  : 


Att'pc  ■■  (l+.v)-  +  (/2//)- 

,,  =  4L  f 

Att'oc  J 


djc 


(l+.r)-'  - 

x-  +  {niff 


As  /}//^  1.  a  good  estimate  is  obtained  from 


(D.lOa) 


(D.lOb) 


<)  f 

-■>  ,  ,  f'O  f  .V  }(}' 

vZ  =  l.  +  L  =  — ; —  d.r ^  =  — (l-7r/4). 

'  -  iTT-pc  I  x-  +  {n/fY  Irr-pc 


(D.ll) 


If  we  define  A„  =  IttcIQ.  the  vacuum  wavelength  of  the  molecular  transition,  and  a  the  lattice  constant, 
we  finally  have  (1  -  7r/4=  1/4), 


{v  I c)-^TT{fin)ia I . 


(D.12) 
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